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ABSTRACT 


Contact  stresses,  contact  forces,  relative  contact 
surface  displacements,  and  the  dissipative  effects  of 
friction  are  computed,  by  the  finite  element  method,  for 
various  two-  and  three-dimensional  contact  problems. 

The  finite  element  technique  is  verified  by  analysis  of 
several  two-dimensional  frictionless  contact  problems; 
the  Hertz  contact  of  two  cylinders,  a  rigid  sleeve/elastic 
shaft  shrinkfit,  and  an  elastic  sleeve/elastic  shaft  shrink- 
fit.  In  these  cases  the  finite  element  calculations  compare 
favorably  to  existing  solutions.  The  contact  analysis 
capability  is  extended  to  frictional  contacts  by  con¬ 
sidering  plane  stress  and  axisymmetric  contact  problems 
with  friction  and  slip,  for  which  reasonable  results  are 
produced.  The  capability  is  further  extended  to  complex 
three-dimensional  contacts  by  an  effort  to  determine  the 
contact  forces  and  frictional  dissipation  taking  place 
in  the  cyclic  bending  of  a  shrinkfit  assembly.  This 
work  demonstrates  that  advanced  nonlinear  finite  element 
methods  can  be  used  to  solve  a  variety  of  mechanical 
engineering  problems  involving  unlubricated  contact 
surfaces  and  the  effects  of  friction. 


ADMINISTRATIVE  INFORMATION 

The  work  described  in  this  report  was  performed  at  the  David  W.  Taylor  Naval 
Ship  Research  and  Development  Center  (DTNSRDC)  under  the  In-house  Research/ In-house 
Engineering  Development  (IR/IED)  program.  Program  Element  61152N,  Task  Area 
ZR0230301,  and  Work  Unit  1720-110. 

METRIC  CONVERSION 

All  numerical  quantities  in  this  report  are  expressed  in  U.S.  customary  units. 
Use  the  following  factors  to  convert  to  metric  units: 

1  in.  ■  2.54  cm 

1  lb  -  0.454  kg 

1  psi  »  0.690  N/cm^ 

1  Ib/in  ■  1.751  N/cm 

1  in-lb  -  0.113  J 
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INTRODUCTION 


In  many  cases,  the  reliability  of  machine  elements  in  power  transmission  and 
power  generation  systems  depends  on  the  integrity  of  contact  surfaces  between  system 
components.  Although  the  wear  resistance  of  component  surfaces  is  a  vital  consid¬ 
eration,  the  performance  of  contact  surfaces  in  a  structural  mechanical  sense  is 
also  very  important.  All  machinery  components  deform  under  load,  and  such  deforma¬ 
tions  influence  both  the  contact  pressures  and  the  frictional  forces  transmitted 
between  contact  surfaces.  These  forces  are  also  influenced  by  the  materials  and 
lubricants  involved. 

Many  examples  can  be  cited  in  which  contact  stresses  and  the  frictional 
behavior  of  contact  surfaces  affect  the  various  possible  failure  modes  of  naval 
ship  and  submarine  machine  components.  Often,  the  contact  surface  is  lubricated 
(e.g.  the  concentrated  contacts  arising  in  gear  and  bearing  assemblies).  Surface 
fatigue  failure  mechanisms  such  as  pitting  and  fretting,  which  occur  in  reduction 
gear  and  propulsion  shaft  bearings,  respectively,  are  influenced  by  contact  pres¬ 
sure  and  frictional  stress.  Contact  pressure  and  friction  resulting  from  component 
deformation  influence  the  performance  of  mechanical  fastenings  (bolts,  joints, 
rivets,  threaded  fasteners)  and  seal  systems  (O-rings  and  gaskets).  The  Influence 
of  contact  stresses  on  failures  is  an  important  consideration  in  high-speed  machin¬ 
ery  components,  such  as  turbine  blade/rotor  fastenings.  Finally,  contact  surface 
gapping  is  sometimes  a  potential  problem  (e.g.  the  behavior  of  preloaded  bolts  and 
blade  palm/hub  surfaces  in  controllable  pitch  propeller  systems). 

When  dealing  with  contact  surface  integrity  problems,  the  navy  engineering 
community  has  traditionally  devoted  much  attention  to  the  materials  aspects. 
Solutions  to  wear,  galling  and  surface  fatigue  problems  have  typically  been  sought 
through  extensive  and  expensive  experimentation  and  materials  evaluations.  The  on¬ 
going  concern  with  submarine  shaft  seals  provides  a  prime  example.  These  machinery 
component  problems  have  been  addressed  almost  exclusively  from  a  materials  stand¬ 
point,  with  relatively  little  attention  to  the  structural  mechanical  performance. 
Clearly,  a  seal  cannot  perform  as  intended  unless  the  component  deformations  and 
the  geometry  of  seal  and  surrounding  components  assure  an  optimum  contact  surface 
load  transmission  pattern.  The  experimental  line  of  attack  was  the  only  available 
recourse  when  analytical  methods  for  contact  and  friction  analysis  were  limited, 
impractical,  or  Impossible.  The  rapid  development  of  finite  element  methods  to 
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analyze  structural  mechanical  systems  of  complex  geometry  and  treat  previously 
Intractable  nonlinear  phenomena  has  changed  this  situation  completely. 

As  naval  machinery  systems  have  become  more  complex,  the  need  for  a  contact/ 
friction  structural  analysis  capability  haB  become  more  acute.  Although  they  retain 
their  importance,  experiments  and  materials  evaluations  that  account  for  the  influ¬ 
ence  of  all  variables  are  difficult  and  expensive  to  conduct.  An  analytical  tool 
is  needed  which  is  capable  of  mathematically  evaluating  the  performance  of  contact 
surfaces  in  mechanical  component  designs.  Such  an  analytical  tool  may  also  be  of 
assistance  in  failure  assessment. 


OBJECTIVE  AND  SCOPE 

The  purpose  of  this  study  is  twofold: 

(1)  To  verify  evaluate,  and  refine  newly  developed  nonlinear  finite 
element  techniques  for  contact  stress  and  friction  analysis. 

(2)  To  apply  these  new  techniques  to  a  complex  practical  engineering 
problem,  fretting  corrosion  fatigue,  in  an  attempt  to  correlate 
calculated  contact  and  friction  behavior  to  experimental  evidence 
of  fretting. 

The  scope  of  this  study  is  limited  in  several  ways: 

(1)  The  theory  is  limited  to  dry,  unlubricated  contacts.  The  important 
problem  area  of  elastohydrodynamic  lubrication  is  not  considered. 

(2)  Of  principal  concern  are  interference  contacts  over  a  large  surface  area; 
concentrated  contacts  are  not  considered. 

(3)  Mixed  contact/gapping  is  not  considered,  although  the  theory  is  capable 
of  straightforward  treatment  of  such  problems. 

(4)  This  study  is  limited  to  static  and  quasi-static  cyclic  loadings,  although 
theories  exist  which  can  handle  dynamic  impact  problems. 

(5)  The  theory  describes  the  effects  of  contact  behavior  only  in  a  macroscopic 
structural  mechanical  sense.  No  attempt  is  made  to  address  the  tribolog¬ 
ical  issues  of  contact  and  friction  behavior  (asperity  interaction,  con¬ 
tact  surface  irregularity,  local  cold  welding  and  cracking,  etc.) 

These  limitations  can  be  placed  in  proper  perspective  by  realizing  that  this 
study  is  only  a  first  attempt  to  apply  state-of-the-art  structural  mechanics  analy¬ 
sis  methods  to  a  wide  class  of  complex,  practical  contact  problems. 
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BACKGROUND 


CLASSICAL  MATHEMATICAL  APPROACHES  TO 
CONTACT  PROBLEMS  OF  SOLID  MECHANICS 

The  mathematical  techniques  of  elasticity  theory  have  long  been  used  to  solve 

the  static  contact  problems  of  solid  mechanics.  The  standard  work  of  Timoshenko 

(reference  1,*  pp.  409-420)  contains  a  treatment  of  the  so-called  Hertzian  contact 

of  spheres  and  cylinders,  the  simplest  type  of  contact  problem  and  the  first  one 

reduced  by  rational  analysis.  The  Hertz  problems  are  also  discussed  by  Lubkln 

(reference  2,  Chapter  42)  and  by  Love  (reference  3,  pp.  192-200).  Many  Russian 

mathematicians  have  specialized  in  contact  research,  and  three  treatises  on  such 

problems  have  been  produced. The  English-language  literature  also  contains  a 

vast  number  of  papers  on  a  wide  variety  of  contact  problems;  a  few  recent  samples 
7-12 

are  referenced. 

The  elasticity  approach  is  important  for  fundamental  understanding  and  has 
provided  many  useful  solutions  to  engineering  problems.  From  the  engineer's  stand¬ 
point,  however,  the  purely  mathematical  approach  presents  three  difficulties: 

(1)  Complexity 

The  commonly  used  specialized  techniques  (integral  transforms,  potential 
functions)  require  a  in-depth  understanding  of  the  theory  of  partial  differen¬ 
tial  equations  and  integral  equations.  As  a  result,  solutions  tend  to  be 
involved  mathematical  expressions  that  are  not  always  presentable  in  closed 
form. 

(2)  Geometrical  limitations 

In  order  for  elasticity  solutions  to  be  tractable,  the  problem  must  be 
posed  in  simple  coordinate  systems  (Cartesian,  radial).  Many  practical 
contact  problems  are  characterized  by  complex  geometry  that  does  not  fit  a 
standard  coordinate  system.  Also,  because  specialized  elasticity  solutions 
contain  assumptions  on  boundary  conditions  and  symmetry  planes,  their  valid¬ 
ity  is  restricted. 

(3)  Limited  threatment  of  nonlinear  effects  in  mathematical  solutions 

When  friction  is  considered, for  example,  (which  is  rare)  the  ideal  assump¬ 
tion  of  full  adherence  is  usually  invoked,  while  the  real  situation  may  involve 
mixed  slip  and  adherence.  In  contact  and  gapping  situations,  the  extent  of 
gapping  must  be  assumed  "a  priori,"  an  assumption  which  may  be  difficult  to 

*A  complete  listing  of  references  is  given  on  page  117. 
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justify.  The  nonlinearity  of  dynamic  contact-impact  kinetics  is  also  very 
difficult  to  treat  analytically. 

This  brief  critique  of  classical  methods  demonstrates  the  desirability  of  an 
approximate  numerical  approach  to  analysis  of  contact  effects  in  complex  naval 
engineering  problems. 

FINITE  ELEMENT  APPROACHES  IN  CONTACT,  FRICTION, 

AND  GAPPING  PROBLEMS  IN  SOLID  MECHANICS 

The  finite  element  method  has  recently  been  extended  to  handle  nonlinear  con¬ 
tact  problems  in  solid  mechanics.  For  many  years,  the  thrust  of  research  has  been 
to  devise  specialized  techniques  for  specific  classes  of  contact  problems.  Earliest 
efforts  addressed  frictionless  contact  arising  from  a  prescribed  indentation  without 
gapping.  J  Techniques  were  then  developed  to  handle  frictionless  contact/gap  prob¬ 
lems  (reference  14,  and  work  by  Gifford*);  frictionless  Interference  fits  without 

gapping;15-*^  and  combinations  of  contact,  gapping,  and  interference.**  The  first 

18-20 

finite  element  treatment  for  contact  and  gapping  with  friction  considered  two- 

dimensional  contacts  modelled  by  low-order  finite  elements.  This  capability  has  been 
extended  to  three-dimensional  problems.  The  general  approach  of  this  work  is  to 
formulate  the  contact  problem  in  terms  of  conventional  incremental  equilibrium  equa¬ 
tions  subject  to  nonlinear  constraint  conditions  on  contact  surface  displacements 
and/or  forces.  This  approach  can  be  conveniently  generalized  in  terms  of  the  clas- 
sicial  Lagrange  multiplier  method.  In  which  nonlinear  contact  constraint  equations 
and  conventional  incremental  equilibrium  equations  are  solved  simultaneously  in  an 
iterative  manner.  The  Lagrange  multiplier  approach  Is  considered  in  this  study. 

Another  line  of  attack  is  to  model  contact  interfaces  as  a  fictitious  layer 
of  material  possessing  empirical  constitutive  properties  that  approximately  describe 
observed  contact  and  friction  behavior.  Contact  and  gap  phenomena  can  be  modelled 
with  "bilinear  springs"  that  possess  very  high  stiffness  when  contact  is  detected 
and  are  assigned  vanishingly  small  stiffness  when  gapping  occurs.  Analogous 
"friction  springs"  can  similiarly  be  activated  and  deactivated  according  to  whether 


Unpublished  work  by  L.N.  Gifford  of  the  Structures  Department  of  DTNSRDC. 

Unpublished  work  by  R.  A.  LIndeman  of  the  U.S.  Naval  Weapons  Laboratory, 
Dahlgren,  Virginia. 
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adherence  or  slip  is  occurring.  *  Contact  surfaces  have  also  been  modelled  as 
a  continuous  layer  of  material  with  bilinear  stiffness. ^  The  bilinear  spring 
approach  has  been  used  in  general-purpose  nonlinear  finite  element  programs. 

These  works  share  a  common  feature;  a  contact  surface  stiffness  parameter  is 
assigned  which  supplements  the  conventional  stiffnesss  coefficients  of  the  incre¬ 
mental  equilibrium  equations.  The  "stiffness"  (also  called  "penalty")  approach 
has  advantages  in  some  situations  and  is  also  considered  in  this  study. 

More  advanced  lines  of  finite  element  contact  research  have  also  been  pursued. 
An  advanced  stiffness-like  method  which  models  complex  frictional  behavior  by  "bond 
elements"  has  been  proposed.  Complex  adherence  and  slip  rules  closely  resembling 
classical  plasticity  theories  have  been  formulated  for  two-  and  three-dimensional 
problems.  »  These  works  model  frictional  behavior  by  rational  constitutive  rela¬ 
tions.  The  constitutive  equation  approach  holds  promise  for  modeling  finer  details 
of  friction  behavior  for  a  wide  variety  of  materials;  however,  it  is  not  convenient¬ 
ly  applicable  to  the  problems  considered  here.  Finally,  a  methoo  has  been  devised 

o  o  oo 

for  finite  element  solution  of  dynamic  contact- impact  problems,  '  which  are 
vitally  important  but  of  little  relevance  to  this  study. 

These  recent  works  have  made  possible  the  approximate  numerical  analysis  of 
a  wide  variety  of  contact  problems.  The  limitations  of  classical  mathematical 
methods,  geometry-dependence,  and  need  for  special  boundary  conditions  are  elimi¬ 
nated  quite  naturally  by  finite  element  techniques.  Highly  nonlinear  contact 
phenomena,  such  as  mixed  contact  and  gapping,  mixed  adherence  and  slip,  and  contact/ 
impact,  are  accomodated  quite  generally  by  finite  element  codes.  Finally,  most  of 
those  engineers  who  do  not  have  the  extensive  mathematical  background  needed  to 
understand  elasticity  methods  can  work  with  even  nonlinear  finite  element  methods. 

The  finite  element  techniques  which  have  been  applied  to  problems  considered 
in  this  study  are  more  fully  described  in  the  following  section. 
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APPROACH 


LAGRANGE  MULTIPLIER  METHOD  FOR  FINITE 
ELEMENT  MODELLING  OF  CONTACT,  GAP, 

AND  FRICTION  PROBLEMS 

The  elastic  stiffness  matrix  and  the  incremental  load  vector  in  the  structural 
incremental  equilibrium  equations  are  formulated  in  this  approach  by  the  generalized 
principle  of  minimum  potential  energy  (or  alternatively,  the  principle  of  virtual 
work).  The  incremental  equilibrium  equations  are  supplemented,  however,  by  con¬ 
straint  terms  of  the  Lagrange  multiplier  type  which  impose  the  contact  or  gap 
conditions  existing  at  the  current  load  step.  Finite  element  meshes  representing 
the  solid  bodies  in  the  problem  are  modelled  in  the  usual  way  by  arrays  of  ele¬ 
ments  interconnected  at  nodal  points.  The  nodal  arrangement  on  surfaces  that 
initially  interfere  (as  in  a  shrinkfit),  initially  contact,  and  initially  contact 
but  gap  apart  as  loading  proceeds,  must  be  such  that  opposing  nodes  are  paired; 
that  is,  their  initial  coordinates  are  either  coincident  (interference  and  contact), 
or  spatially  close  together  and  in  line  with  one  another  (initial  gapping,  see 
Figure  1).  It  is  assumed  that  displacements  in  the  problem  are  small,  so  that 
opposing  nodes  which  are  initially  close  together  remain  so  in  the  course  of  load¬ 
ing.  The  analyst  specifies  "contact  elements"  or  "contact  kinematic  constraints" 
between  chosen  pairs  of  candidate  contact  nodes.  The  relative  normal  displacements 
are  monitored  at  these  node  pairs  as  the  structure  is  loaded.  The  Lagrange  multi¬ 
pliers,  which  are  calculated  unknowns  as  are  the  nodal  displacements,  turn  out  to 
be  the  normal  contact  forces  transmitted  between  the  node  pairs.  The  initial 
interference  case,  which  exists  if  all  node  pairs  are  overlapping,  is  a  linear 
problem  which  can  be  solved  in  one  step  since  all  Lagrange  multipliers  are  active. 

In  the  general  case,  however,  the  contact  conditions  change  in  the  course  of 
loading  (i.e.  gaps  close  or  open)  and  the  Lagrange  multipliers  are  made  active 
(gap  closed,  contact  force  transmitted)  or  set  to  zero  (gap  open,  no  contact  force), 
depending  on  the  current  normal  relative  displacements.  This  constitutes  the  non¬ 
linearity  of  the  problem.  A  complete  mathematical  description  of  the  Lagrange 
multiplier  theory  and  the  associated  finite  element  equations  are  given,  in  vec- 
torial  form,  by  Hibbitt  and  Rubin.  Specialized  examples  of  finite  element  equa¬ 

tion  systems  have  been  developed  from  the  general  forms  in  reference  22.  These 


unpublished  developments*  Illustrate  the  adherence/s  lip  constraints  and  the  solution 
process  for  various  contact  conditions  of  practical  interest* 

The  Lagrange  multiplier  approach  is  also  used  to  model  the  frictional  behavior 

22 

of  contact  surfaces*  Classical  Coulomb  friction  theory  is  applied;  this  assumes 
that  an  adherence  condition  exists  at  a  contact  node  pair  when  the  local  frictional 
force  (force  tangential  to  the  contact  surface)  is  below  the  local  static  frictional 
force  limit.  Slip  (relative  tangential  displacement)  occurs  when  the  frictional 
force  reaches  this  limit.  The  friction  limit  is  defined  by  an  empirical  relation 
of  the  normal  contact  force  Fjj,  in  which  U  is  the  experimentally  measured  static 
friction  coefficient: 

Fy  <  p  •  Fjj  (adherence)  (1) 

FT  “  m  '  fN  (slip)  (2) 

In  the  complete  adherence  condition,  all  relative  tangential  displacements  on 
the  contact  surface  are  constrained  to  zero,  and  the  frictional  forces  are  unknown. 
In  the  mixed  adherence  and  slip  condition,  the  frictional  force  is  constrained  to 
Fj  -  U  •  Fn  at  slipping  node  pairs,  and  the  slip  is  the  unknown  quantity  associated 
with  those  node  pairs.  The  displacement  response  to  external  loading  is  nonlinear 
in  a  load  step;  this  causes  the  Coulomb  slip  limit  to  be  reached  at  some  node  pairs 
while  the  others  remain  adhered.  The  nonlinear  problem  is  solved  by  monitoring 
the  relation  of  normal  force  to  frictional  force,  the  appropriate  constraints  being 
imposed  according  to  whether  equations  (1)  or  (2)  hold.  After  a  trial  solution 
that  assumes  full  adherence  everywhere,  final  equilibrium  at  the  current  load  is 
found  through  an  iterative  process  that  unlocks  adhered  nodes  not  satisfying  the 
Coulomb  limit. 

The  Lagrange  multiplier  approach  for  both  normal  contact/gapping  and  friction 

behavior  has  been  incorporated  as  the  "contact  elements"  of  the  MARC  general  purpose 
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nonlinear  finite  element  code.  Several  finite  element  contact  approaches 

are  very  similar  to  the  MARC  method  but  do  not  contain  a  frictional  slip  degree  of 

freedom.  They  can  be  considered  to  be  special  cases  of  the  MARC  approach.  As  it 

turns  out,  the  Lagrange  multiplier  method  works  quite  well  for  the  normal  contact 

* 

Authored  by  D.  E.  Lesar  of  the  Structures  Department  of  DTNSRDC* 


force/gap  part  of  the  contact  problem  but  does  not  always  work  well  for  friction. 
Hibbitt,  the  developer  of  this  method,  has  Indicated  that  the  Lagrange  multiplier 
approach  for  friction  and  slip  suffers  mathematical  instabilities  and  leads  to  con¬ 
vergence  difficulties  in  cases  where  adherence  to  slip  transitions  occur. ^  The 
method  appears  to  predict  adherence/slip  behavior  reasonably  well  in  the  2-D  analy- 
sis,  as  long  as  the  extent  of  slip  is  small.  Attempts  to  apply  MARC  to  the  3-D 
friction  problems  considered  here  resulted  in  severe  convergence  problems.  This 
poor  performance  led  to  the  consideration  of  a  new  friction  theory  which  borrows 
from  the  so-called  "penalty  method"  used  in  finite  element  solutions  of  contact 
problems.  This  improved  friction  theory  is  discussed  in  the  next  section. 

STIFFNESS  METHOD  FOR  FINITE  ELEMENT  MODELLING 
OF  FRICTIONAL  BEHAVIOR  IN  CONTACT/GAP 
AND  FRICTION  PROBLEMS 

The  penalty  (or  stiffness)  method  has  been  used  in  the  past  to  model  the 
contact/gap  aspects  of  contact  problems.  Candidate  contact  node  pairs  are  linked 
by  what  amount  to  simple  truss  elements  that  possess  a  "bilinear  stiffness."  These 
elements  are  assigned  a  very  high  stiffness  when  the  current  nodal  displacements 
indicate  a  closed  gap  and  a  very  low  stiffness  when  gapping  occurs.  This  method, 
used  in  ANSYS^  and  a  previous  version  of  MARC,^  tends  to  predict  physically 
invalid  contact  surface  overlap  when  the  "stiffnesses”  are  not  optimally  tuned 
to  the  particular  problem  at  hand.  This  poor  performance  led  to  development  of 
the  Lagrange  multiplier  approach  discussed  earlier. 

Some  penalty  method  concepts  have  been  resurrected  in  development  of  a  new 

friction-constraint  element  that  circumvents  possible  convergence  problems  often 

encountered  in  the  Lagrange  multiplier  friction  approach.  This  feature,  contained 
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in  the  ABAQUS  general-purpose  nonlinear  finite  element  program  ,  retains  Lagrange 
multiplier  contact/gap  modelling  but  makes  use  of  a  completely  new  friction  approach 
described  briefly  in  the  following  paragraphs. 

The  ABAQUS  theory  assumes  that  the  frictional  shear  force  magnitude  FT  is 
linearly  proportional  to  the  tangential  contact  surface  displacement  by  a  "fric¬ 
tional  stiffness"  or  "stiffness  in  stick”  parameter  as  long  as  the  frictional 
shear  force  is  less  that  the  static  frictional  limit.  Because  of  K^,  the  contacting 
solids  behave  as  if  connected  by  elastic  springs  that  are  capable  or  transmitting 
frictional  force.  These  "springs"  represent,  in  a  gross  way,  the  resistance  of  the 
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contact  surface  to  shear  force  and  can  be  Interpreted  as  the  Integrated  effect 
of  the  microscopic  asperity  interaction  occurring  in  reality.  This  approach,  which 
allows  relative  tangential  dispacements  on  the  contact  surface,  is  more  realistic 
than  the  Lagrange  multiplier  method,  which  assumes  that  all  such  displacements 
are  zero  in  the  adherence  regime.  A  danger  exists,  however;  must  be  carefully 
chosen  so  as  to  produce  neither  unrealistic  results  nor  numerical  instability  prob¬ 
lems.  This  matter  is  discussed  in  the  analysis  section. 

The  slip  condition  is  defined  in  ABAQUS  by  a  frictional  shear  force  limit  similar 

to  equation  2.  Instead  of  calculating  Lagrange  multipliers  associated  with  slip, 

ABAQUS  directly  imposes  the  frictional  slip  constraint  on  the  frictional  forces  by 
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a  so-called  "radial  return"  algorithm. 

ABAQUS  regards  an  incremental  solution  as  "converged"  when  all  nonequillbrated 
(residual)  nodal  forces  resulting  from  the  approximate  nature  of  the  piecewise  linear 
solution  to  the  nonlinear  problem  are  within  some  selected  tolerance.  This  conver¬ 
gence  criterion  is  quite  stringent  and  is  particularly  amenable  to  reliable  solution 
of  contact  and  friction  problems  where  the  nonlinearities  are  related  to  forces. 

Although  finite  element  equations  for  the  ABAQUS  contact  and  friction  approach 
are  completely  developed  in  vectorial  form,  some  simple  specific  examples  written 
in  matrix  form  provide  additional  insight.  Lesar  of  the  DTNSRDC  Structures  Depart¬ 
ment  has  prepared  such  developments  for  simple  examples  that  parallel  two  cases 
described  for  Lagrange  multiplier  friction  in  unpublished  notes  mentioned  previously. 

CONTACT  CONSTRAINT/CONTINUUM  ELEMENT  COMPATIBILITY 

Both  the  Lagrange  multiplier  contact  "elements"  of  MARC  and  the  hybrid  Lagrange 
multiplier-penalty  method  contact  "elements"  of  ABAQUS  are  kinematic  constraints 
imposed  at  specified  contact  surface  node  pairs.  Both  are  inherently  compatible 
with  all  displacement-based  continuum  solid  finite  elements  of  low  interpolation 
order.  Such  pointwise  constraints  are,  however.  Incompatible  with  second-order 
solids.  Additional  constraints  are  needed  to  enforce  compatibility. 

The  compatibility  difficulty  can  be  clarified  by  considering  the  simple  solid 
element/contact  element  combinations  shown  in  Figures  2a  and  2b.  In  a  finite  ele¬ 
ment  assembly,  distributed  pressure  loads  are  transmitted  from  one  solid  element  to 
another  via  equivalent  nodal  forces,  which  are  calculated  according  to  the  element's 
displacement  interpolation  functions.  Low-order  elements  (4-node  axisymmetric, 
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8-node  3-D  solid)  will  convert  a  uniform  pressure  load  imposed  on  one  face  to  a  set 
of  equal  compressive  nodal  forces  on  the  opposite  face.  In  Figure  2a,  the  contact 
elements  transfer  these  compressive  forces  to  the  rigid  foundation.  Second-order 
elements  (8-node  axisymmetric ,  20-node  3-D  solids)  convert  this  same  uniform  pres¬ 
sure  load  to  a  set  of  both  compressive  and  tensile  nodal  forces  on  the  opposite 
face  (Figure  2b,  see  also  reference  39,  p.  223).  As  a  result,  the  contact  elements 
at  outer  corners  are  forced  to  carry  tensile  load,  an  obviously  unreasonable  result 
that  leads  to  false  predictions  of  gapping  at  these  corner  nodes.  Similar  diffi¬ 
culties  arise  with  tangential  displacement  degrees  of  freedom,  resulting  in  invalid 
frictional  force  and  slip  predictions.  The  advantages  offered  by  higher-order  solid 
elements  in  modelling  complex  geometries  are  lost  unless  this  compatibility  problem 
is  remedied. 

The  simplest  cure  is  to  convert  the  quadratic  displacement  interpolation  func¬ 
tions  for  contact  surface  variables  to  linear  functions.  The  equivalent  nodal 
forces  generated  by  these  linear  functions  will  always  have  the  same  algebraic  sign. 
Linearization  is  accomplished  by  condensing  all  contact  surface  midside  nodes  out  of 
the  equilibrium  equation  system,  i.e.,  by  imposing  constraints  of  the  form  shown  in 
Figure  3: 


u£  =  1/2  (0*  ~  Up  (3A) 

=  1/2  (U*  -  Up  (3B) 

for  axi  ymmetric  elements  and 

u£  =  1/2  (U®  ~  Up  (4A) 

u£  -  1/2  (Uj  -  U^)  (4B) 

u£  =  1/2  (U®  "  Up  (4C) 


for  three  dimensional  solids. 

The  MARC  and  ABAQUS  programs  contain  features  allowing  easy  input  of  such 
constraints.  However,  the  local  linearization  remedy  possesses  two  undesirable 
features.  First,  part  of  the  advantage  of  high-order  elements  is  lost.  The 
removal  of  midside  nodes  reduces  the  number  of  contact  node  points.  Loss  of  con¬ 
tact  surface  variables  degrades  the  attainable  accuracy.  Second,  linearization  is 
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strictly  valid  only  when  the  contact  surface  is  initially  composed  of  straight 
lines.  In  cases  of  curved  geometry,  the  displacement  linearization  destroys  the 
smoothness  of  predicted  deflection  distributions  (Figure  4).  This  error  worsens 
with  greater  element  arc  length  and  greater  contact  surface  curvature. 

The  linearization  constraint  was  the  only  compatibility  remedy  available  when 
the  analyses  reported  herein  were  carried  out.  Compatibility  difficulties  have  been 
solved,  however,  through  development  of  a  new  concept  in  contact  and  friction  ele¬ 
ments  for  the  ABAQUS  program.  The  basis  for  these  new  elements  is  a  reformulation 
of  contact  and  friction  constraints  into  a  continuously  interpolated  form,  in  much 
the  same  manner  in  which  solid  element  properties  are  approximated.  Fami1ies  of 
first-  and  second-order  "interface  elements"  for  two-  and  three-dimensional  analy¬ 
sis  have  been  devised  which  interpolate  contact  and  frictional  forces  with  chosen 
functions  that  are  inherently  compatible  with  first-  or  second-order  two-  and  three- 
dimensional  solids.  The  new  ABAQUS  approach requires  use  of  solid  elements 
with  variable  numbers  of  nodes.  The  extra  data  preparation  price  for  these  special 
elements  is  well  worth  the  much  improved  accuracy  and  generality  which  results. 

ANALYSIS 

VERIFICATION  OF  LAGRANGE  MULTIPLIER  METHOD 
FOR  STATIC  FRICTIONLESS  CONTACT 

The  normal  contact  force,  gapping,  and  friction  capability  of  the  Lagrange 
multiplier  method  has  been  partially  verified  in  previous  work.  The  MARC  program 
was  used  to  analyze  the  Hertzian  contact  of  two  infinitely  long  cylinders  of  equal 
radius.  The  predicted  contact  pressure  and  internal  stress  distributions  agreed 
reasonably  well  with  the  classical  solution  (reference  1,  pp.  414-420).  An  inter¬ 
ference  fit  contact  of  two  cylinders  with  a  friction  coefficient  of  .30  was  also 

analyzed.  Frictional  slip  at  all  nodes  resulted  from  imposition  of  a  sufficient 

22 

horizontal  rigid  body  displacement  on  one  cylinder. 

An  elasticity  solution  of  a  similar  fully  slipped  contact  problem^  results  in 
an  internal  stress  distribution  that  is  nonsymmetric  with  respect  to  the  center 
plane  of  the  contact  zone.  MARC  predicted  stress  contours  which  were  very  similar 
to  the  elasticity  solution. 22 

Additional  comparisons  to  elasticity  solutions  were  thought  necessary,  mainly 
to  clear  up  solid  eleraent/contact  constraint  compatibility  issues  but  also  to  gain 
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a  feel  for  the  behavior  of  contact  constraints  in  progressively  more  complex 
idealizations  of  mechanical  systems*  One  such  system  of  interest  to  the  Navy's 
ship  power  transmission  community  is  the  shrinkfit  of  a  finite  length  sleeve  onto 
a  long  shaft.  The  principal  focus  of  this  study  is  on  shaft  and  sleeve  contact 
problems,  the  first  step  being  evaluation  of  MARC  solutions  for  stresses,  contact 
forces,  and  displacements  in  two  axisymmetric  shrinkfit  problems  which  have  been 
treated  by  elasticity  methods. 

Axisymmetric  Rigid  Sleeve  and  Elastic  Shaft  Shrinkfit 

Case  1.  A  rigid  sleeve  of  finite  length  shrunk  onto  a  long,  hollow,  elastic  shaft 
is  analyzed.  The  contact  is  frictionless  and  the  shrinkfit  is  a  uniform  Indentation 
over  the  entire  contact  zone.  This  elasticity  solution^  is  more  complete,  in 
terms  of  calculated  results,  than  many  shrinkfit  studies.  The  physical  situation 
and  important  dimensions  of  the  problem  are  shown  in  Figure  5.  A  schematic  of 
the  finite  element  representation  appears  in  Figure  6.  Second-order  (8-node)  axi¬ 
symmetric  solids  were  used  with  linearization  of  contact  surface  displacements. 

The  shrinkfit  interference  is  specified  as  a  negative  initial  gap.  A  geometric 
symmetry  plane  through  the  center  of  the  contact  zone  is  used  to  advantage.  The 
boundary  conditions  of  the  elasticity  solution  can  only  be  approximated.  The  in¬ 
finitely  long  shaft  cannot  be  modelled  without  special  "exponential  decay”  elements 
(reference  39,  pp.  660-664)  so  a  finite  shaft  half-length  eight  times  the  contact 
zone  half-length  was  chosen.  The  end  of  the  shaft  is  assumed  to  be  free  to  expand; 
as  will  be  seen,  this  was  an  inappropriate  assumption.  A  rotation-free  boundary 
cannot  be  directly  specified  at  the  symmetry  plane  since  rotation  degrees  of  freedom 
are  not  available  for  continuum  elements.  The  standard  "roller"  boundary  condition 
was  used  and  this  effectively  eliminated  almost  all  rotation.  The  finite  element 
mesh  is  pictured  in  Figure  7;  the  finer  discretization  at  the  edge  of  the  contact 
zone  (an  area  of  steep  stress  gradients)  is  evident. 

It  should  be  apparent  that  this  problem  could  be  solved  without  aid  of  contact 
elements  -  one  may  well  have  simply  input  the  shrinkfit  as  a  displacement  load. 

This  case  served  to  resolve  some  side  issues  relevant  to  contact  modelling.  One 
main  thrust  of  this  effort  was  to  see  if  MARC  could  produce  reasonable  stress  pre¬ 
dictions  in  a  fairly  simple  shaft  contact  problem.  A  second  objective  was  to  eval¬ 
uate  the  possibility  of  embedding  layers  of  linear  (4-node  axisymmetric)  elements 


within  the  8-node  element  mesh  in  the  contact  region.  This  seemed  to  be  a  viable 
way  of  achieving  contact  constraint/solid  element  compatibility.  A  third  objective 
was  to  prove  that  contact  surface  linearizations  are  really  needed  in  a  practical 
problem. 

The  third  issue  was  settled  quite  easily  by  an  analysis  without  linearization, 
which  produced  a  scalloped  indentation  pattern  instead  of  the  required  straight 
line;  an  obvious  consequence  of  the  algebraic  sign  problem  with  equivalent  nodal 
forces.  A  second  analysis  used  a  pure  8-node  element  mesh  and  a  third  analysis  con¬ 
tained  4  layers  of  4-node  elements  as  shown  in  Figure  8.  Element  mixing  turned  out 
to  be  a  poor  option  because  linear  elements  cannot  capture  steep  stress  gradients; 
a  typical  result  being  the  badly  discontinuous  integration  point  stresses  shown  in 
figure  9.  The  mixed  mesh  solution  also  resulted  in  higher  computer  charges.  The 
element  mixing  option  seemed  to  be  a  poor  performer  and  was  discarded. 

The  correlation  of  MARC  results  to  the  findings  of  Hill  et  al.^  is  judged  to 
be  fair  despite  two  difficulties.  First,  Hill's  numerical  results  are  for  an  in¬ 
compressible  material  with  a  Poisson's  ratio  of  1/2.  This  number  can  be  used 
validly  in  finite  element  calculation  only  if  special  incompressible  elements  are 
used.  Conventional  elastic  elements  with  V  »  .49  were  used  here.  Second,  the  effect 
of  infinite  length  may  have  been  better  represented  with  a  fixed  boundary  at  the 
end  of  the  shaft;  a  boundary  free  of  restraint  was  used  instead.  This  choice  may 
have  produced  some  disagreement  in  radial  (Figures  10a,  10b)  and  axial  displacements 
(Figure  11).  The  predicted  trends  are,  however,  correct. 

Computed  variation  of  the  four  stress  components  with  shaft  length  is  compared, 
for  several  depths  below  the  contact  surface,  to  the  Hill  et  al.  solution  in  Fig¬ 
ures  12-15.  The  stress  curves  are  "by  eye"  interpretations  of  integration  point 
values;  no  attempt  was  made  to  extrapolate  and/or  smooth  them  mathematically.  The 
correlation  appears  fairly  good,  despite  the  Poisson's  ratio  approximation.  Pre¬ 
liminary  analyses  with  meshes  coarser  than  the  one  pictured  in  Figure  8  showed 
severe  overshoot  and  oscillation  in  radial,  axial,  and  hoop  stress  beyong  the  edge 
of  the  rigid  sleeve.  This  effect  was  minimized  by  the  finer  mesh  used  to  obtain 
results  shown  in  Figures  10-15. 

It  must  be  emphasized  that  the  elasticity  solution  predicts  singular  radial 
and  hoop  stress  at  the  edge  of  the  contact  zone,  a  typical  feature  of  contact 
problems  containing  corner  discontinuities.  Finite  element  models  require  use  of 
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special  displacement  functions  (reference  39,  pp.  664-673)  when  singular  stress 
behavior  needs  to  be  accurately  simulated.  The  sharp,  finite  peaks  predicted  by 
MARC  are  not  physically  real  and  there  is  little  hope  of  predicting  the  actual 
stress  state  in  the  vicinity  of  a  contact  corner  unless  inelastic  deformations 
are  allowed.  At  this  point,  it  is  clear  that  a  truly  realistic  finite  element 
solution  (in  terms  of  stresses)  for  a  contact  problem  with  singularities  Involves 
many  complicated  side  issues  which  were  pursued  to  some  extent  in  the  next  analysis. 

An  error  in  modelling  was  made  but  not  identified  until  much  later,  in  that 
contact  elements  were  specified  at  both  corner  and  midside  node  pairs.  This  should 
have  resulted  in  a  computational  error  since  radial  displacements  at  midside  nodes 
are  constrained  out  of  the  problem.  Despite  this,  MARC  reported  compressive  forces 
at  all  node  pairs.  A  contact  pressure  distribution  can  be  computed  by  dividing 
these  forces  by  areas  tributary  to  the  node  pairs.  This  can  be  done  for  all  node 
pairs  (both  corner  and  midside)  individually  or  by  allocating  midside  node  pair 
contact  forces  to  corner  nodes  and  using  tributary  areas  that  are  twice  as  large 
as  in  the  first  case.  Both  options  predict  roughly  the  same  pressure  distribution. 
The  only  significant  difference  is  that  the  individual  node  pair  option  leads  to 
an  edge  pressure  that  is  twice  as  large  as  the  midside  node  pair  allocation  option. 
The  pressure  distribution  shown  in  Figure  16  corresponds  to  the  second  (allocation) 
option.  This  is  believed  to  be  the  proper  way  of  interpreting  the  computed  contact 
forces  under  the  erroneous  modelling  conditions.  Note  that  these  pressures  agree 
fairly  closely  with  the  radial  stresses  at  integration  points  just  below  the  surface 
(see  Figure  12) . 

In  summary,  the  MARC  program  proved  itself  capable  of  modelling  the  essential 
features  of  a  fairly  simple  shrinkfit  problem,  although  the  predicted  stress  pat¬ 
terns  are  of  questionable  accuracy  due  to  boundary  condition  errors,  Poisson’s 
ratio  effects,  and  influence  of  a  contact  pressure  singularity. 

Axisymmetric  Elastic  Sleeve  and  Elastic  Shaft  Shrinkfit 

Case  2.  An  elastic  sleeve  of  finite  length  shrunk  onto  a  long,  solid  elastic  shaft 
is  analyzed.  The  contact  is  frictionless  and  the  shrinkfit,  as  specified,  is  uni¬ 
form;  however,  in  reality  it  is  nonuniform  because  of  the  elasticity  of  both  members. 
The  physical  situation  is  shown  in  Figure  17,  and  a  schematic  of  the  finite  element 
representation  with  dimensions  appears  in  Figure  18.  The  modelling  scheme,  as  far 
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as  boundary  conditions,  shrinkfit  loading,  etc.  are  concerned,  is  exactly  as  in 
Case  1  except  that  dimensions  are  different  and  the  sleeve  must  be  modelled.  The 
material  is  taken  to  be  steel  (E  -  30x10^  psi,  u  *.3). 

This  problem  is  a  stronger  test  of  contact  element  performance.  The  dimensions 
chosen  correspond  to  a  shrinkfit  assembly  subjected  to  fretting  tests  at  DTNSRDC.* 
This  particular  assembly  was  one  of  a  series  subjected  to  cyclic  bending  tests  in 
an  effort  to  identify  key  structural,  material,  and  loading  variables  leading  to 
fretting  fatigue  failure.  (The  relation  of  contact  analysis  to  fretting  is  post¬ 
poned  to  later  discussion;  the  primary  purpose  here  is  to  discuss  contact  modelling 
performance.)  This  geometry  does  not  correspond  to  any  set  of  proportions  con- 
sidered  by  Conway  and  Farnhani0  but  the  MARC  results  can  be  bracketed  by  two  elas¬ 
ticity  solutions.  A  more  comprehensive  mathematical  shrinkfit  study^  produced 
results  which  correspond  only  to  shorter  sleeves  but  are  helpful  for  qualitative 
comparison  purposes. 

The  finite  element  mesh  used  here  is  shown  in  Figure  19.  It  became  apparent  in 
Case  1  that  a  mesh  convergence  study  in  conjunction  with  a  smoothing/extrapolation 
scheme  is  necessary  to  get  the  best  attainable  estimate  of  stresses  very  close  to  or 
at  the  contact  surface,  subject  to  limitations  implied  by  the  inherent  singularity. 
The  mesh  convergence  issue,  although  important,  seemed  less  worthy  of  attention 
than  the  stress  smoothing  question. 

It  is  well  known  that  the  finite  element  procedure  based  on  assumed  displace¬ 
ment  fields  gives  continuous  and  smooth  displacement  variations  but  yields  stresses 
that  are  discontinuous  between  elements.  This  mathematical  fact  of  life  has  been 
dealt  with  through  many  semi-empirical  schemes  (nodal  averaging,  extrapolation) 
for  finding  smoothed  stress  fields.  Rational  means  of  smoothing  have,  however, 
been  devised;^  ^  one  such  method  is  used  here.  In  second-order  elements  (8-node 
axi8ymmetric,  20-node  solids)  each  stress  component  local  to  each  element,  when 
integrated  exactly,  follows  a  parabolic  form.  These  parabolic  forms  are,  under 
some  situations,  fairly  continuous  from  element  to  element  but  tend  to  become  badly 
behaved  as  stress  gradients  steepen  (Figure  20).  Use  of  "reduced  integration" 


These  tests  were  conducted  by  W.  Werchniak  of  the  Ship  Materials  Engineering 
Department  of  DTNSRDC.  The  results  of  the  tests,  herinafter  called  "DTNSRDC  tests," 
have  not  been  published. 
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elements  (In  which  variables  are  Integrated  Inexactly  with  fewer  Integration  points) 
Is  equivalent  to  a  least-squares  smoothing  of  the  local  element  stress  field. ^ 

The  smoothed  local  stresses  follow  a  linear  variation  between  integration  points, 
and  smoothed  stresses  can  be  found  anywhere  on  an  element  edge  by  extrapolation. 

A  fairly  realistic  estimate  of  nodal  point  stress  can  then  be  found  after  additional 
extrapolation  concluded  by  averaging  (Figure  21).  This  method  cannot  be  claimed  as 
the  "best,"  but  it  proved  both  simple  and  expedient  and  produced  reasonable  results 
in  this  study. 

The  MARC  displacement  predictions  are  compared,  in  Figures  22-23,  to  solutions^ 

that  bracket  the  finite  element  results  in  terms  of  nondimenslonal  contact  zone 

size.  Both  the  radial  displacement  and  the  percentage  of  the  specified  shrinkflt 

prevented  by  the  shaft  fall  between  the  elasticity  solutions  and  follow  similar 

patterns.  The  nonuniform  radial  contact  surface  deflection  predicted  by  MARC  corre- 

« 

sponds  fairly  well  to  the  uniform  plane  strain  prediction  from  the  classical  Lame 

/  P 

shrinkflt  solution.  The  shaft  prevents  85-90  percent  of  the  unrestrained  sleeve 
shrinkage.  The  finite  element  results  would  fall  halfway  between  the  bracketing 
curves  if  the  sleeve  length-to-thickness  parameter  b/h  corresponded  to  the  b/h 
considered  by  Conway  and  Farnham.  The  slightly  lower  MARC  b/h  lowers  the  percentage 
shrinkflt  restraint  somewhat.  The  shaft  restraint  is  greater  at  the  end  of  the 
sleeve,  so  this  effect  is  less  pronounced  there.  Note  that  even  though  local  dis¬ 
placement  linearizations  are  in  effect  on  the  contact  surface,  the  nodal  displace¬ 
ments  could  still  be  faired  into  the  smooth  curves  of  Figures  22-23.  The  axial 
displacements  on  the  shaft  surface  and  inner  sleeve  bore  are  shown  in  Figure  24; 
this  agrees  qualitatively  with  previously  cited  results. The  sleeve  bore  moves 
in  opposition  to  the  shaft,  the  shaft  expanding  outward  and  the  sleeve  contracting 
inward. 

The  reduced  integration/extrapolation  smoothing  method  discussed  earlier  was 
used  to  obtain  estimates  of  radial,  hoop,  and  axial  stress  on  the  contact  surface. 

The  results  of  this  effort  are  shown  in  Figures  25-27.  Radial  stress  estimates 
are  computed  from  integration  points  in  the  sleeve  as  well  as  integration  points 
in  the  shaft.  The  estimates  based  on  shaft  data  agree  almost  perfectly  with  the 
plane  strain  stress  prediction^®  and  also  fall  between  existing  bracketing  cases. ^ 
The  sleeve  extrapolation  is  in  lesser  agreement  because  the  stress  gradients  are 
much  steeper  there,  and  the  mesh  is  not  fine  enough  to  model  this  gradient 
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accurately.  Two-thirds  of  the  contact  zone  Is  essentially  In  plane  strain  while 
Interesting  edge  effects  In  the  form  of  "wiggles”  are  evident  (the  Appendix  contains 
a  brief  discussion  of  why  these  wiggles  exist).  The  singular  stress  at  the  sleeve 
edge  Is  evident.  Although  the  singularity  appears  much  less  severe  than  in  the 
rigid  sleeve/elastic  shaft  case,  the  stress  peaks  at  the  sleeve  edge  cannot  be 
construed  as  realistic  predictions. 

The  axial  stress  behavior  agrees  qualitatively  with  other  results. ^  The 
sleeve  load  causes  shaft  expansion,  which  is  restrained  somewhat  by  the  long 
shaft  beyond  the  sleeve  edge.  The  resulting  compression  is  counteracted  by 
sleeve  tension  and  shaft  tension  outside  the  contact  zone.  This  effect  is  highly 
local  to  the  edge,  however.  Finally,  the  shear  stress  extrapolated  to  the  con¬ 
tact  surface  is  essentially  zero,  as  it  should  be  for  a  frictionless  contact. 

An  evaluation  of  contact  element  performance  is  shown  in  Figure  28.  All 
contact  elements  reported  a  closed  condition  and  bore  compressive  forces.  Since 
this  ideaMzation  also  had  contact  elements  at  constrained  midside  node  pairs,  con¬ 
tact  pressures  were  calculated  from  contact  forces  by  allocating  midside  node  pair 
forces  to  adjacent  corner  node  pairs.  As  in  Case  1,  this  does  not  seem  to  have 
seriously  affected  the  results.  The  extrapolated  radial  stresses  should  approach 
these  calculated  contact  element  pressures  in  the  limit.  Figure  28  provides 
proof  that  contact  element  predictions  are  in  good  agreement  with  extrapolated 
stress,  except  at  the  edge  of  the  contact  zone.  Here,  contact  elements  show  a 
tendency  toward  a  sharp  stress  peak  while  extrapolated  stresses  do  not.  This  is 
a  fault  of  the  solid  elements  rather  than  contact  elements.  The  solids  do  not 
have  the  requisite  enhancement  for  modelling  singular  stress  fields. 

In  summary,  the  finite  element  method  of  contact  modelling  through  Lagrange 
multipliers  worked  admirably  well  in  two  axisymmetric  shrinkfit  problems.  Intei^- 
element  compatibility  can  be  attained  by  linearizing  contact  surface  displacements. 
Although  many  difficult  side  issues  arise  when  accurate  contact  surface  stresses  are 
desired,  the  problem  can  be  solved  If  enough  effort  is  devoted  to  mesh  convergence 
Issues,  and  smoothing  and  extrapolation  of  predicted  stresses.  Calculated  stresses 
near  singularities  must  be  viewed  with  great  caution  in  such  problems. 
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EVALUATION  OF  LAGRANGE  MULTIPLIER  METHOD 
FOR  3-D  STATIC  FRICTIONLESS  CONTACT 

The  methodology  of  contact  modelling  is  now  extended  to  three  dimensions  by 
reanalysis  of  the  axisymmetric  elastic  shaft/elastic  sleeve  problem.  The  primary 
purposes  here  are  to  verify  the  linearization  compatibility  enforcement  method  for 
higher  order  3-D  (20-node)  elements  and  to  evaluate  the  effect  of  a  coarser  mesh. 

The  modelling  scheme,  as  far  as  boundary  conditions,  loads,  and  linearizations  are 
concerned, is  much  like  Case  2.  The  shrinkfit  radial  interference  is  .001  in.  at 
all  contact  node  pairs,  and  the  material  is  again  assumed  to  be  steel.  Reduced 
integration  elements  are  again  used,  and  contact  suface  stresses  are  obtained  by 
the  same  extrapolation  method.  Two  symmetry  planes  are  utilized,  as  shown  in  Figure 
29,  with  the  X  *  0,  Y  =*  0,  and  2  **  0  planes  fixed  against  X,  Y,  and  Z  displacement 
respectively.  This  mesh  is  very  similar  to  that  used  in  an  earlier,  independent 
attempt  to  analyze  the  DTNSRDC  fretting  fatigue  test  rig  with  the  NASTRAN  program.* 
The  original  intent  was  to  study  the  response  of  this  contact  idealization  to  the 
.0005  in.  radial  shrinkfit  and  nonsymmetric  loads  considered  in  the  NASTRAN  analy¬ 
sis,  but  the  MARC  program  appeared  to  be  prohibitively  expensive  for  this.  The 
NASTRAN  work  did  not  consider  the  frictionless  axisymmetric  response  so  the  analy¬ 
ses  could  not  be  compared,  anyway. 

The  calculated  results  agree  quite  closely  with  those  of  the  axisymmetric 
model  considered  previously.  Radial  displacements  predicted  by  the  axisymmetric 
and  3-D  idealizations  are  compared  in  Figure  30.  The  linearizations  did  not  cause 

O 

any  significant  circumferential  bias;  that  Is,  radial  displacements  at  6  «  45  are 

O 

hardly  different  from  those  at  8  “  90  .  The  details  of  edge  effect  on  stresses 
and  contact  element  pressures  differed  because  of  the  coarser  mesh  near  the  outer 
part  of  the  contact  zone.  This  effect  is  seen  In  Figure  31,  in  which  smoothed  and 
extrapolated  nodal  averages  of  radial  stress  are  compared  for  the  axisymmetric  and 
3-D  models.  The  stress  curves  shown  have  been  faired  through  nodal  averages  based 


The  nonlinear  effects  of  contact  and  frictional  slip  were  treated  by  succes¬ 
sive  linear  analyses  interrupted  by  analyst  Intervention.  Contact  node  pairs  were 
first  constrained  to  move  together.  If  contact  surface  stresses  showed  that  separa¬ 
tion  or  slip  should  be  occurring,  appropriate  constraints  were  released  and  the 
analysis  redone  until  gaps  were  eliminated  and  shear  stresses  were  below  or  close 
to  the  frictional  limit.  This  unpublished  work  was  conducted  by  E.  Schroeder  of 
the  Computation,  Mathematics  and  Logistics  Department  of  DTNSRDC. 
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on  shaft  Integration  point  data.  Although  the  stress  decay  beyond  the  edge  of 
the  contact  zone  is  much  slower  for  the  coarser  3-D  mesh,  the  essentials  of  the 
edge  effect  are  captured  and  plane  strain  conditions  are  correctly  approached. 

Radial  stresses  based  on  sleeve  extrapolation  are  different,  but  it  appears  that 
better  agreement  could  be  obtained  from  sleeve  mesh  refinement;  this  is  evident 
from  Figure  32.  Since  contact  elements  were  erroneously  specified  once  again  at 
constrained  midside  node  pairs,  a  special  scheme  allocating  these  forces  to  corner 
nodes  had  to  be  devised  (Figure  33).  This  adjustment  resulted  in  a  reasonable 
contact  pressure  distribution  that  agreed  quite  well  with  the  pressure  prediction 
of  the  axisymmetric  idealization  (Figure  34).  Not  surprisingly,  the  edge  pressures 
differed  substantially  due  to  the  different  discretization  levels.  The  mistake  of 
retaining  contact  elements  at  midside  nodes  was  corrected  in  later  analyses. 

In  summary,  the  3-D  analysis  results  agreed  with  the  axisymmetric  idealization 
results,  despite  obvious  discretization  deficiencies.  The  contact  elements  per¬ 
formed  acceptably  well  in  a  3-D  model  with  higher-order  elements  when  contact 
surface  displacement  linearizations  were  used. 

EVALUATION  OF  LAGRANGE  MULTIPLIER  METHOD  FOR 
STATIC  CONTACTS  WITH  FRICTION  AND  SLIP 

The  Lagrange  multiplier  method  of  MARC  produced  only  limited  success  in  model¬ 
ling  frictional  contact  behavior.  Two  attempts  worked  reasonably  well;  the  axisym¬ 
metric  elastic  shaf t/elastic  sleeve  shrinkfit  with  friction  included,  and  a  very 
coarsely  discretized  3-D  model  of  an  asymmetrically  loaded  shrinkfit.  MARC  failed 
to  produce  convergent  solutions  to  a  plane  stress  contact,  friction,  and  slip  prob¬ 
lem  and  a  more  elaborate  model  of  the  3-D  asymmetric  shrinkfit.  These  four  cases 
are  discussed  in  the  following  section. 

Axisymmetric  Elastic  Sleeve  and  Elastic  Shaft  Shrinkfit 

Case  1 .  The  elastic  sleeve/elastic  shaft  model  is  reanalyzed  with  a  friction 
coefficient  of  0.15.  The  problem  is  completely  linear  because  contact  surface 
shear  forces  are  not  large  enough  to  cause  slip.  As  a  consequence,  all  contact 
node  pairs  remain  in  adherence.  The  only  variables  affected  to  any  discernible 
extent  by  friction  were  the  contact  surface  axial  displacements,  the  axial  stresses, 
and  the  shear  stresses.  In  Figures  35  and  36,  results  for  the  first  two  variables 


are  compared  to  the  frictionless  predictions.  The  relative  slip,  or  difference 
between  shaft  and  sleeve  movement,  is  less  than  in  the  frictionless  case,  a  physi¬ 
cally  reasonable  result.  Friction  also  converts  axial  stress  from  a  localized 
edge  effect  to  an  effect  over  almost  the  entire  plane  strain  zone,  a  natural  conse¬ 
quence  of  frictional  shear  constraint  on  the  contact  interface.  The  smoothed  inte¬ 
gration  point  shear  stresses  were  so  drastically  discontinuous  from  element  to 
element  that  various  efforts  to  further  smooth  the  integration  point  data  could 
not  produce  consistent  results.  Apparently,  the  shear  stress  gradient  in  the  radial 
direction  is  so  steep  that  only  contact  element  shear  forces  are  numerically  signi¬ 
ficant.  A  frictional  shear  stress  distribution  was  calculated  from  contact  element 
shear  forces  in  the  same  manner  as  for  normal  contact  pressures,  by  reallocating 
midside  node  pairs  (Figure  37).  Although  the  trend  toward  zero  shear  at  the  symmetry 
plane  and  a  peak  at  the  sleeve  edge  is  very  rough,  Figure  37  looks  somewhat  like 
the  fully  adhered  frictional  shear  predictions  for  short  sleeve/long  shaft  shrink- 
fits.^  The  lesson  learned  is  that  symmetry  boundary  contact  node  pairs  should 
be  assigned  zero  friction  coefficient  to  suppress  spurious  frictional  shear.  This 
correction  was  made  in  subsequent  analyses.  In  this  case,  the  frictional  shears 
are  very  small  compared  to  contact  pressures;  their  ratio  is  at  most  of  order 
1/100,  much  less  than  the  specified  Coulomb  slip  limit  of  0.15.  Interestingly,  the 
result  of  Conway  and  Farnham's^  plane  strain  punch/slab  frictional  contact  study 
shows  that  the  extent  of  slip  in  such  a  shrinkfit  depends  only  on  geometric  ratios 
and  not  on  material  properties.  It  appears  that  slip  would  occur  in  this  problem 
if  the  sleeve  were  much  shorter. 

2-D  Friction  and  Slip  Problem 

Case  2 .  The  previous  case  was  particularly  simple  since  no  frictional  slip  occurred. 
In  a  more  thorough  test  of  Lagrange  multiplier  friction  theory,  approximate  adher¬ 
ence  and  slip  zones  for  a  rigid  flat-ended  punch  indenting  a  planar  elastic  slab 
are  calculated  by  a  combined  elasticity/numerical  method. ^  Conway  and  Farnham's 
calculation  shewed  that  the  adherence/sllp  zone  size  depends  only  on  the  Coulomb 
friction  coefficient  and  the  punch  half-width  to  slab  depth  ratio.  Although  the 
load  in  Figure  2  in  Conways ’s  work  appears  as  a  concentrated  force,  the  unknown 
contact  surface  pressures  and  shears  are  determined  by  assuming  a  uniform  indenta¬ 
tion  over  the  whole  contact  zone,  a  condition  forced  by  the  rigid  punch.  The  MARC 


solutions  were  carried  out  with  a  Young's  Modulus  of  30  x  10^,  Poisson's  ratio 
of  zero  and  a  uniform  contact  surface  Indentation  of  .0005  In.  The  infinitely 
long  slab  is  approximated  by  a  finite  slab  with  stress-free  ends  that  are  free 
to  displace;  the  affect  of  a  fixed  end  was  found  to  be  negligibly  small.  The 
finite  element  mesh  is  shown  in  Figure  38;  eight-node  reduced  integration  plane 
stress  elements  are  used  with  the  usual  contact  surface  displacement  linearizations. 
The  Young's  Modulus  and  indentation  were  varied  in  some  test  runs  but  these  vari¬ 
ables  had  no  affect  on  adherence/slip  predictions.  The  MARC  calculations  covered 
the  friction  coefficient  range  .10  <_  y  <_  .60  for  both  a/h  =>  1/4  and  a/h  =  1/2. 

The  solution  strategy  was  the  same  as  for  all  previous  analyses: 

(i)  Impose  shrinkage  in  increment  ’’zero,"  assume  zero  friction,  resolve 
contact/gapping  by  iteration,  find  normal  contact  forces. 

(ii)  Impose  friction  in  increment  “one,"  update  normal  contact  forces 
retaining  same  shrinkage  or  indentation,  resolve  adherence/slip 
conditions,  find  tangential  contact  (friction)  forces. 

(iii)  Repeat  (ii)  until  MARC  convergence  criteria  are  satisfied. 

The  MARC  results  are  summarized  in  Table  1.  The  only  valid  results  agreeing 
with  reference  49  were  for  complete  adherence  at  the  largest  friction  coefficients. 
None  of  the  other  solutions  are  valid;  MARC  either  repeated  a  nonconvergent  solution 
or  ended  up  In  a  slowly  convergent  iteration  loop  which,  for  the  lowest  friction 
coefficients,  approached  a  totally  invalid  result  (full  adherence  where  nearly  com¬ 
plete  slip  Is  the  correct  answer).  Efforts  to  locate  a  conceptual  error  or  a 
modelling  blunder  were  fruitless.  This  failed  attempt  led  to  further  investigations 
of  the  MARC  friction  capability. 

Primitive  Model  of  3-D  Shaft  and  Sleeve  Shrinkfit 
with  Monotonic  Bending 

Case  3.  The  MARC  friction  capability  is  now  tested  In  a  very  coarsely  discretized 
shaft/sleeve  shrinkfit  subjected  to  nonsymmetric  bending  load.  The  shrinkfit 
assembly  is  the  same  as  that  treated  in  the  elastic  shaft/elastic  sleeve  problem 
considered  earlier.  The  finite  element  idealization  is  shown,  with  boundary  condi¬ 
tions,  in  Figure  39.  Note  that  only  two  symmetry  planes  can  be  utilized.  The 
radial  shrinkfit  interference  is  .005  in.  and  the  bending  load  magnitiude  is  static¬ 
ally  increased  from  zero  to  1000  in-lb.  A  friction  coefficient  y  -  .15  is  specified. 
The  dimensions  and  load  parameters  pertained  to  a  shrinkfit  assembly  in  DTNSRDC 


tests  that  suffered  fretting  corrosion  damage  when  subjected  to  a  cyclic  1000  in-lb 
bending  load.  This  assembly  is  pictured  in  Figure  40.  (The  relation  of  contact 
analysis  to  the  fretting  experiments  is  deferred  to  a  later  section.  It  suffices 
to  say  at  this  point  that  adherence/slip  predictions  are  important.) 

Nodal  forces  equivalent  to  the  bending  load  are  calculated  in  a  preliminary 
MARC  analysis  of  the  shaft  alone.  The  axial  deflections  corresponding  to  a  1000 
in-lb  moment  are  calculated  from  a  linear  elastic  beam  equation  analysis  of 
the  shaft.  A  MARC  shaft  model  is  fixed  at  one  end  with  these  axial  displacements 
imposed  at  the  free  end.  The  axial  reaction  forces  reported  by  MARC  at  the  fixed 
end  are  equivalent  to  the  desired  moment.  The  maximum  lateral  deflection  obtained 
is  very  close  to  the  simple  beam  theory  prediction.  The  sign  of  the  moment  is 
immaterial  since  the  shaft  is  geometrically  symmetric  about  the  neutral  plane  of 
bending. 

After  bending  moment  nodal  forces  were  identified,  MARC  solutions  of  shaft/ 
sleeve  interaction  were  carried  out.  A  series  of  MARC  analyses  intended  to  find  the 
load  incrementation  scheme  necessary  for  convergence  showed  that  nine  equal  bending 
moment  increments  were  needed.  Convergence  could  not  be  obtained  for  three  or  six 
increments.  The  finite  element  model  is  so  coarse  that  calculated  stresses  have 
little  meaning,  but  deflections  and  the  frictional  behavior  of  the  contact  surface 
are  of  interest.  MARC  predicted  a  maximum  laterial  shaft  deflection  of  about  .00221 
In. ,  which  agreed  closely  with  .00224  in.  predicted  by  an  approximate  stepped  beam 
deflection  formula. Deflections  in  the  plane  of  bending  showed  that  the  assembly 
bends  as  a  classical  thick  beam  with  noticable  shear  deflection  in  the  sleeve. 

The  slip  and  adherence  history  is  shown  schematically  in  Figure  41.  The 
shrinkfit  Is  slightly  relieved  on  one  side  and  slightly  increased  on  the  diametri¬ 
cally  opposed  side  in  the  plane  of  bending,  a  reasonable  result.  No  gapping  occurs 
for  any  load;  the  bending  moment  is  not  great  enough  to  relieve  the  initial  shrink- 
fit  Interference.  The  ratios  of  frictional  force  to  normal  force  and  the  percentage 
of  shrinkfit  restrained  by  the  shaft  which  are  attained  at  each  contact  surface 
node  pair  are  summarized  In  Table  2  for  the  equilibrated  shrinkfit  condition  and 
various  levels  of  bending  moment.  Frictional  slip  begins  soon  after  bending  load 
is  applied,  and  the  zone  of  slip  spreads  somewhat  with  further  application  of 
load.  Slip  is  confined  to  the  upper  and  lower  edges  of  the  contact  zone  in  the 
plane  of  bending.  This  behavior  is  not  unexpected  and  the  MARC  friction  theory 
seems  to  have  worked  well  in  this  particular  case. 
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Improved  Model  of  3-D  Shaft  and  Sleeve  Shrinkfit 
with  Monotonic  Bending 

Case  4.  The  same  modelling  techniques  used  successfully  in  Case  3  are  now  applied 
to  the  same  shaft/sleeve  interaction  problem  with  Improved  discretization  in  the 
radial  and  axial  directions.  The  new  finite  element  mesh  is  shown  in  Figure  42. 

All  contact  surface  displacements  are  linearized,  and  contact  elements  link  only 
the  corner  node  pairs. 

The  idealization  is  well-behaved  for  zero  friction  and  solution  for  normal 
contact  force  for  the  bending  moment  range  0-1000  in-lb  was  easily  obtained.  The 
friction  analyses  failed,  however,  to  approach  a  reasonable  result.  All  attempts 
to  find  load  increments  small  enough  for  convergence  failed  due  to  slow  convergence, 
divergence,  or,  more  commonly,  "ping-pong"  divergence?  These  problems  arose  when 
slip  initiated  at  one  or  more  contact  surface  node  pairs.  Load  increments  as 
small  as  1/216  of  maximum  load  (4.63  in-lb)  managed  to  isolate  the  slip  of  one 
node  pair  at  a  time  but  did  not  eliminate  "ping-ponging,”  which  occurred  for  steps 
of  1/9,  1/18,  1/36  and  1/72  of  maximum  load  as  well.  Apparently,  at  initial  slip 
the  first  estimates  of  friction-force  to  current  normal  force  ratio  are  so  much 
greater  than  the  coefficient  of  friction  that  convergence  is  impossible  no  matter 
how  small  the  load  increment.  In  one  interesting  analysis,  for  example,  one  node 
pair  surpasses  the  Coulomb  limit  during  a  load  increment  of  4.63  in-lb  (1/216  maxi¬ 
mum  load).  MARC  wound  up  in  a  "ping-pong”  loop  with  the  tangential-to-normal  force 
ratio  changing  from  .103  to  .547  and  back  to  .103,  etc.,  never  coming  close  enough 
to  the  actual  limiting  value  of  0.15. 

A  thorough  check  of  Input  by  plotting  and  line-by-line  scanning  of  input  data 
failed  to  reveal  errors  or  inconsistencies.  Consultation  with  MARC  users  and  devel¬ 
opers  did  not  reveal  any  faults  in  the  problem  definition.  The  unsatisfactory 
performance  of  MARC  in  this  case  and  in  the  2-D  punch  problem  prompted  consideration 
of  an  alternative  computer  code  for  frictional  adherence/slip  behavior.  At  this 
time,  the  ABAQUS  code  became  available  for  in-house  use.  It  appeared  to  be  worth¬ 
while  to  try  the  ABAQUS  program,  which  contains  a  contact  and  friction  capability. 

^'Ping-Pong"  divergence  can  occur  when  an  iterative  nonlinear  equation  solution 
process  fails  to  home  in  on  a  unique  equilbrium  state  at  the  end  of  a  load  increment. 
In  the  particular  cases  encountered  here,  the  iterations  bounced  from  one  totally 
different  "solution"  to  another,  neither  of  which  fully  satisfied  force  equilibrium 
or  contact  constraints. 


The  contact  algorithms  of  ABAQUS  make  use  of  the  same  Lagrange  multiplier  approach 
used  successfully  in  MARC,  and  the  friction  agorithms  are  based  on  a  stiffness 
me  thod . 

VERIFICATION  AND  EVALUATION  OF  FRICTIONAL 
STIFFNESS  METHOD  FOR  STATIC  FRICTIONAL 
ADHERENCE  AND  SLIP  PROBLEMS 

The  ABAQUS  program  was  applied  to  the  2-D  mixed  adherence/slip  problem  that 
MARC  failed  to  handle  and  to  the  coarsely  discretized  shaft/sleeve  shrinkfit,  bend¬ 
ing,  and  friction  problem.  It  produced  satisfactory  solutions  in  both  cases,  which 
are  described  in  the  following  section. 

2-D  Friction  and  Slip  Problem 

Case  1.  The  mixed  adherence/slip  contact  problem^  is  modelled  by  ABAQUS  in  much 
the  same  way  as  in  MARC.  Only  the  a/h  =  1/2  case  is  considered;  the  mesh  is 
identical  to  that  seen  in  Figure  38.  The  specified  indentation,  material  proper¬ 
ties,  and  dimensions  are  the  same  as  used  in  MARC;  the  same  elements  (8-node 
reduced  integration  plane  stress)  are  used  along  with  the  same  contact  surface 
displacement  linearizations.  The  one  new  item  of  input  is  the  "stiffness  in  stick" 
or  "frictional  stiffness,"  a  parameter  discussed  earlier  in  the  section  on  stiffness 
approaches  to  contact.  This  can  be  interpreted  as  the  ratio  of  the  local  contact 
surface  shear  force  to  the  local  contact  surface  tangential  displacement.  However, 
because  this  parameter  is  problem  dependent,  the  rationale  for  its  estimation  is 
not  at  first  obvious.  The  friction  theory's  developer  says  that  the  proper  stiff¬ 
ness  to  use  is  the  highest  number  for  which  a  convergent  solution  is  possible. 

The  Invalid  MARC  results  for  this  problem  provided  ratios  of  friction  force  to 
tangential  contact  displacement  and,  with  no  other  guidance  available,  sufficed 
as  a  basis  for  a  parametric  study.  For  u  =  0.35  these  ratios  fell  within  the  range 
3  x  10^  <  Kj.  <  5  x  10^  when  all  contact  node  pairs  are  considered.  The  results 
of  a  successful  ABAQUS  parametric  study  based  on  these  numbers  is  shown  in  Table 
3  for  U  =  0.35  and  Table  4  for  y  =  0.80.  The  upper  limit  on  frictional  stiffness 
is  very  sharply  defined  by  a  convergence  limit  while  the  lowest  stiffness  results 
in  a  total  adherence  condition.  Intermediate  values  produce  mixed  adherence/slip 
solutions.  Total  adherence  Is  not  reasonable  for  y  =  0.35,  hence  the  lower  can 
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be  discarded.  For  U  =  0.80,  total  adherence  appears  to  be  practically  the  only 
possible  solution,  which  is  somewhat  consistent  wii_h  previously  cited  findings. ^ 
Additional  mesh  refinement  and  trials  are  needed  to  obtain  the  optimum  Kf,  but 
this  is  hardly  worth  the  effort.  Note  that  for  U  -  0.80,  the  convergence  limit  is 
again  sharply  defined. 

From  the  above  two  studies,  it  appears  that  -  1.1  x  10^  is  adequate;  further 

ABAQUS  simulations  were  carried  out  with  this  chosen  stiffness  parameter  for  other 

friction  coefficients.  Resulting  ratios  of  contact  forces  are  shown  in  Table  5,  and 

49 

the  ABAQUS  adherence  zone  predictions  are  compared  to  previously  cited  research 
in  Figure  43.  Although  the  results  are  rather  approximate  because  of  discretiza¬ 
tion  limitations  and  the  inexactness  of  Kf,  the  correct  trends  are  predicted  quite 
adequately. 

The  stiffness  approach  to  friction  modelling  requires  a  frictional  stiffness  K^, 
which  may  be  estimated  if  some  guidance  is  already  available.  The  difficulty 
in  the  stiffness  method  is  to  determine  a  physically  reasonable  K.£  in  a  completely 
unique  problem.  This  issue  is  addressed  in  the  next  case. 

Primitive  Model  of  3-D  Shaft  and  Sleeve 
Shrinkfit  with  Monotonic  Bending 

Case  2.  The  ABAQUS  friction  capability  is  now  applied  to  the  coarsely  discretized 
3-D  shaft/sleeve  problem  shown  in  Figure  39.  The  mesh,  element  type,  material 
properties,  and  loading  (.005  in.  unrestrained  sleeve  shrinkage  followed  by  1000 
in-lb  bending  load)  are  the  same  as  before.  All  contact  surface  displacements  are 
restricted  to  linear  variations.  The  friction  coefficient  is  fixed  at  U  =  0.15, 
and  frictional  stiffness  is  varied  in  a  parametric  study  to  find  reasonable 
limits  for  this  particular  problem. 

A  preliminary  analysis  of  frictionless  axisymmetric  contact  without  bending 
produced  normal  contact  forces  within  3.5%  of  the  MARC  predictions.  The  ABAQUS 
model's  maximum  lateral  displacement  was  almost  identical  to  the  MARC  value.  Subse¬ 
quent  analysis  focused  on  the  frictional  stiffness  issue. 

The  behavior  of  this  system  with  respect  to  differing  is  much  more  complex 
than  in  the  2-D  rigid  punch/elastic  slab  problem.  All  contact  nodes  in  the  punch 
problem  bear  increasing  normal  forces  as  the  indentation  increases.  In  the  3-D 
problem,  however,  some  contact  node  pairs  carry  increasing  normal  force  while  other 
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node  pairs  are  relieved  of  load  as  bending  proceeds.  Friction  forces  tend  to 
increase  faster  than  normal  forces  at  points  bearing  increasing  normal  load,  and 
normal  forces  decrease  more  slowly  than  friction  forces  at  points  where  normal 
load  is  shed.  As  a  result,  slip  occurs  at  both  locations.  The  friction  response 
to  this  nonproportional  loading  is  clearly  sensitive  to  frictional  stiffness  choice 
in  some  unknown  way. 

A  successful  MARC  solution  to  the  3-D  problem  was  available  as  guidance  for 
the  K£  choice.  "Frictional  stiffnesses"  were  taken  to  be  the  ratio  of  friction 
force  to  relative  tangential  displacement  from  the  MARC  results  for  the  frictional 
shrinkfit  with  no  bending.  All  such  computed  "stiffnesses"  fell  in  the  range 
1  x  10^  <  Kj  <  5  x  10^.  A  series  of  ABAQUS  analyses  for  this  range  of  frictional 
stiffness  at  every  contact  node  pair  was  compared  to  the  MARC  results  for  all 
loads  up  to  the  maximum  bending  load.  A  partial  comparison  of  frictional/normal 
force  ratios  in  Table  6  shows  that  the  solutions  compare  well  initially  but  diverge 
as  bending  load  increases.  Additional  ABAQUS  runs  at  higher  K£  showed  that  the 
force  ratios  agreed  with  MARC  results  only  at  higher  Kf's,  and  the  K^'s  necessary 
for  good  correspondence  increased  with  the  bending  load.  An  example  for  one  node 
pair  is  shown  in  Figure  44.  The  MARC  approach  effectively  determines  the  changing 
contact  surface  frictional  stiffnesses  as  load  increases;  this  stiffness  apparently 
changes  over  orders  of  magnitude  in  a  highly  nonlinear  manner.  The  ABAQUS  approach 
enforces  a  constant  input  stiffness.  If  the  guess  is  too  low,  as  in  this  case, 
frictional  slip  will  not  be  predicted.  If  the  stiffness  is  too  high,  slip  will 
occur  at  too  small  a  relative  displacement.  The  results  shown  in  Table  6  occurred 
because  the  initial  MARC  stiffnesses  are  too  small  to  predict  the  slip  behavior 
that  occurs  as  the  assembly  is  bent. 

Several  preliminary  ABAQUS  analyses  for  order(s)  of  magnitude  variations  in 
were  needed  over  the  entire  range  of  loads.  The  results  of  a  lengthy  parametric 
study  are  partially  summarized  in  Figure  45.  ABAQUS  predicts  practically  all  pos¬ 
sible  adherence/slip  combinations  within  the  range  1  x  10^  <  <  1  x  10^.  Again, 

the  lower  limit  is  a  complete  adherence  solution.  The  higher  extreme  is  certainly 
a  convergence  limit,  as  found  in  the  2-D  punch/slab  problem,  but  this  particular 
limit  was  not  determined. 

The  friction  force  response  of  selected  contact  node  pairs  for  various  choices 
of  frictional  stiffness  is  shown,  for  adherence  conditions,  in  Figure  46  and  for 
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mixed  adherence/s  lip  In  Figure  47.  The  simultaneous  loading  and  unloading  of  con¬ 
tact  node  pairs  and  the  effect  of  slip  on  subsequent  frictional  response  Is  clearly 
seen.  The  Interaction  of  friction  and  contact  forces  Is  indeed  an  interesting  and 
complex  phenomenon. 

The  load  versus  calculated  relative  tangential  displacement  at  an  adhering  con¬ 
tact  node  pair  and  a  slipped  node  pair  for  =*  2  x  10^  is  shown  in  Figure  48.  The 
nonlinear  response  to  monotonic  load  is  easily  explained  if  the  correspondence  of 
the  stiffness  friction  formulation  to  the  White-Besseling  model  of  plasticity^ » ^ 
is  recognized.  The  stiffness  theory  characterizes  the  frictional  response  of  an 
individual  contact  node  pair  by  two  fixed  parameters;  a  "stiffness  in  stick,”  or 
frictional  stiffness,  and  a  coefficient  of  friction.  The  friction  coefficient  de¬ 
fines  a  normal  force-dependent  limit  on  frictional  force  at  each  node  pair.  The 
frictional  stiffness  is  analagous  to  the  elastic  modulus  and  the  limiting  frictional 
force  corresponds  to  the  yield  stress  in  the  elastic-perfectly  plastic  model  of 
metal  plasticity.  In  the  White-Besseling  concept,  a  plastic  material  "element" 
exhibiting  strain-hardening  is  characterized  by  a  suitably  chosen  collection  of 
elastic-perfectly  plastic  "subelements."  The  gross  "element,"  composed  of  many  sub¬ 
elements,  will  exhibit  strain-hardening  when  subject  to  monotonically  increasing 
load.  In  the  ABAQUS  friction  theory,  the  contact  surface  is  characterized  by  a 
suitably  chosen  collection  of  frictional  node  pairs,  each  with  a  force-deformation 
law  analogous  to  an  elastic-perfectly-plastic  White-Besseling  "subelement. ’’  The 
nonlinear  hardening  response  seen  in  Figure  48  Is,  therefore,  not  surprising  and 
can  be  Interpreted  to  demonstrate  a  "contained  plastic  flow"  of  sorts,  the  "flow" 
being  confined  to  the  contact  surface.  An  Interesting,  simple  contact  surface 
stiffness  model  explaining  experimentally  observed"^  dissipative  behavior  in  mono¬ 
tonically  and  cyclically  loaded  contacts  is  found  in  a  journal  article  by  Burdekin, 
et  al.^^  Although  the  authors  utilize  nonlinear  normal  and  frictional  stiffnesses, 
they  mistakenly  characterize  the  phenomena  as  being  "elastic." 

The  friction  theory  is  complicated  by  the  fact  that  the  Coulomb  limit  depends 
on  the  current  normal  force.  It  is  clear,  however,  that  the  frictional  response 
at  any  contact  node  pair  depends  on  both  the  normal  force  and  adherence/sllp  condi¬ 
tions  existing  elsewhere  on  the  surface.  The  strong  resemblance  of  this  friction 
behavior  under  monotonic  loading  to  contained  plastic  flow  at  small  strains  prompts 
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the  question  of  whether  ABAQUS  can  predict  dissipative  frictional  behavior  under 
fully  reversed  cyclic  loading.  This  is  the  focus  of  further  analyses. 

The  ABAQUS  contact/friction  model  works  fairly  reliably  and  is  capable  of  pro¬ 
ducing  physically  reasonable  results.  The  calculations  are,  however,  only  as  good 
as  the  chosen  frictional  stiffness.  The  above  parametric  study  identified  a  reason¬ 
able  range  which  can  only  be  justified  by  engineering  judgment  although  one 
simple  theory^-5  implies  that  reflects  some  measure  of  contact  surface  asperity 
roughness.  The  lowest  for  which  one  node  pair  slipped  in  the  entire  range  of 
loading  is  2  x  10-’,  and  the  lowest  Kf  for  which  three  node  pairs  slipped  during  the 
shrinkfit  is  6  x  10^.  The  low  stiffness  is  easily  justifiable  since  it  is  quite 
near  the  full  adherence  limit.  The  high  stiffness  is  more  subjective  since  it  is 
at  least  two  orders  of  magnitude  less  than  the  limiting  for  convergence,  which 
was  never  determined  here.  These  two  limits  were  used  in  a  subsequent  study  of 
this  coarsely  discretized  system's  reponse  to  cyclic  bending  loads. 

APPLICATION  OF  FRICTIONAL  STIFFNESS  METHOD 
TO  A  COMPLEX  3-D  FRETTING  CORROSION 
AND  FRETTING  FATIGUE  PROBLEM 

The  previously  discussed  MARC  analyses  have  verified  the  Lagrange  multiplier 
method  for  contact  problems,  and  the  ABAQUS  analyses  have  verified  the  stiffness 
approach  for  friction  and  slip.  These  contact  analysis  capabilities  have  been 
extended  to  3-D  problems,  and  the  calculation  of  stresses  near  contact  surfaces  has 
been  proven  possible,  with  some  limitations  due  to  stress  singularities.  The  veri¬ 
fied  and  extended  ABAQUS  finite  element  capability  is  now  applied  to  fretting  corro¬ 
sion  and  fretting  fatigue,  a  complex  and  difficult  class  of  mechanical  engineering 
problems. 

Fretting  corrosion  is  a  type  of  surface  damage  that  results  from  small  periodic 
relative  motions  between  metal  parts  that  are  held  together  by  clamping  pressure. 
Such  conditions  exist  in  many  machine  components  that  are  not  intended  to  undergo 
relative  movement,  e.g.  bolted  or  riveted  connections,  and  shrinkfitted  shaft/sleeve 
assemblies.  Fretting  corrosion  can  drastically  reduce  the  fatigue  strength  of 
machine  parts.  If  it  is  severe  enough,  fretting  corrosion  leads  to  surface  and/or 
near-surface  fatigue  crack  formation,  which  can  ultimately  result  in  crack  propaga¬ 
tion  into  the  bulk  material  and  subsequent  failure  by  fracture.  Fretting  corrosion 
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damage  mechanisms,  from  a  metallurgical  and  tribological  point  of  view,  is  fully 
discussed  by  Duquette.-’*’  Both  this  source^  and  Nishioka,  et  al.^  include  lists 
of  references  on  fretting  corrosion  tests  made  over  the  past  40  years. 

The  study  of  tribological  and  material  influences  on  fretting  fatigue  is  very 
important  and  constitutes  the  bulk  of  related  work  to  date.  Experimental  evidence 
shows,  however,  that  fretting  fatigue  in  machine  components  also  depends  on  macro¬ 
scopic  factors,  some  of  which  can  be  quantified  by  continuum  (structural)  analysis. 
The  relative  slip  amplitude  at  the  edge  of  the  contact  surface  in  pressfitted  hub/ 
axle  assemblies  loaded  In  periodic  (rotating)  bending  was  measured  experimentally 
and  reported  upon.^  These  experiments  showed  that  the  degree  of  relative  slip 
depends  on  geometric  characteristics  (axle  diameter  vs.  hub  diameter,  hub  overhang 
vs.  no  overhang).  Relative  slip  increased  with  increasing  shear  force  and  nominal 
bending  stress,  implying  a  dependence  on  external  loads.  Relative  slip  quickly 
assumed  a  steady-state  variation  which  persisted  for  perhaps  1000  cycles  but  then 
slowly  decreased  in  amplitude  as  wear  processes  began.  The  slip  was  found  to  be 
independent  of  the  rapidity  of  cycling  and  the  material's  surface  hardness.  The 
relation  between  nominal  bending  stress  and  relative  slip  took  the  form  of  a  hys¬ 
teresis  loop,  indicating  the  dissipative  nature  of  the  friction  process.  The  dis¬ 
sipated  energy  is  partly  conducted  and/or  radiated  away  from  the  contact  surface 
as  heat.  The  remaining  energy  is  dissipated  in  plastic  deformation  and  crack  forma¬ 
tion  local  to  the  contact  surface. 

The  DTNSR DC  experiments  attempted  to  define  the  effects  of  contact  interfer¬ 
ence,  sleeve  diameter,  bending  load  amplitude,  and  number  of  bending  cycles  on 
fretting  corrosion  damage  and  fretting  fatigue  failure  in  shrunkfit  shaft/sleeve 
assemblies  loaded  In  rotating  bending.  Fretting  damage,  when  it  occurred,  was 
localized  to  the  edge  of  the  shaft/sleeve  contact  zone.  Fretting  corrosion  showed 
very  little  dependence  on  duration  of  test  or  degree  of  shrinkfit  interference. 
Although  fretting  damage  tended  to  increase  somewhat  for  larger-d iameter  sleeves, 
the  extent  of  damage  depended  more  strongly  on  bending  moment  amplitude  than  any 
other  variable;  fretting  damage  increasing  with  greater  bending  moment.  The  most 
highly  loaded  specimens  eventually  failed  by  fracture. 

Both  the  DTNSRDC  tests  and  the  study  of  Nishioka  et  al.^  indicate,  to  varying 
degrees,  that  fretting  damage  and  relative  slip  are  dependent  upon  geometry  and 
external  loads.  These  variables  can  be  easily  accomodated  in  a  continuum  structural 
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analysis,  particularly  a  finite  element  method  that  handles  contact  and  friction. 
Quasistatic  analysis  Is  appropriate  since  relative  slips  appear  to  be  Insensitive 
to  the  rapidity  of  loading,  and  the  extent  of  fretting  damage  depends  only  slightly 
on  the  number  of  load  cycles.  The  previously  observed^  change  of  relative  slip 
with  cycles  is  due  to  wear  effects  which  cannot  yet  be  modelled  by  finite  elements. 

The  dependence  of  fretting  on  loads  and  geometry  Implies  that  the  stresses 

on  the  contact  surface  are  influential  variables.  This  study  has  shown  that  finite 

element  contact  analysis  methods  can  be  used  to  determine  these  stresses.  However, 

sharp  contact  zone  edges  cause  mathematical  singularities  which  complicate  the 

stress  field  considerably.  If  the  location  of  contact  surface  cracking  is  known 

(as  it  is  here)  and  the  degree  of  initial  crack  can  be  characterized,  linear  elastic 

fracture  mechanics  concepts  can  be  used  to  determine  stress  intensity  factors. 

These  factors  can,  in  turn,  be  correlated  to  crack  growth  rate.  Such  an  approach 

has  been  taken  in  a  simple  analysis  where  the  stressing  is  caused  by  both  general 
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structural  loading  and  by  local  shearing  contact  stresses.  A  similar  approach 
is  taken  for  a  cracked  sheet  subjected  to  predefined  normal  and  shearing  "fretting 
forces."  In  any  case  a  suitable  means  of  stress  singularity  representation  is 
required  to  model  the  singular  stress  field  due  to  special  corner  geometries  or 
existing  cracks.  The  difficulties  with  singularities  detracts  somewhat  from  the 
utility  of  correlating  calculated  contact  surface  stresses  to  fretting  fatigue 
failure. 

Additional  analytical  refinements  beyond  the  scope  of  this  study  are  required 
to  predict  fretting  fatigue  failure;  however,  the  contact  analysis  method  may  be 
useful  in  predicting  fretting  corrosion  damage.  The  finite  element  method  is 
capable  of  predicting  both  the  location  and  extent  of  frictional  slip  in  a  mono- 
tonically  loaded  static  contact.  The  method  should  also  be  capable  of  predicting 
the  tangential  contact  force  vs.  slip  response  in  cyclic  loading.  If  such  dissi¬ 
pative  behavior  (the  referenced  hysteresis  loops'^)  can  be  predicted,  then  the 
energy  converted  into  plastic  deformation  and  heat  generation  local  to  the  contact 
surface  can  be  calculated.  The  damaging  energy  is  that  which  cannot  be  conducted 
and/or  radiated  away  as  heat.  Since  experimental  evidence  shows  that  the  fretting 
process  is  time-independent  (not  accounting  for  wear  effects  on  slip  amplitude), 
a  quasistatlc  analysis  seems  feasible.  Two  series  of  shaft/sleeve  analyses  under 
cyclic  bending  conditions  have  been  completed.  It  is  important  to  realize  that 
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cyclic  bending  is  not  the  same  mode  of  loading  as  the  rotating  bending  conditions 
in  the  DTNSRDC  tests;  this  fact  is  clarified  in  Figure  49.  The  finite  element 
predictions  quoted  in  the  following  sections  do  not  correspond  to  the  real  case 
because  load  is  removed  from  all  points  simultaneously  in  cyclic  bending.  A  zero¬ 
loading  condition  does  not  exist  in  the  actual  experimental  mode  of  rotating  bend¬ 
ing.  Rotating  bending  requires  a  full  three  dimensional  analysis,  however,  without 
symmetry  planes  of  any  kind.  This  was  judged  to  involve  too  much  sophistication 
and  expense  at  this  point  since  the  effects  of  mesh  refinement  and  the  validity 
of  an  assumed  frictional  stiffness  with  respect  to  load  level  were  not  sufficiently 
clear.  For  these  reasons,  efforts  were  confined  to  the  simpler  cyclic  bending  case. 
All  ABAQUS  analyses  discussed  in  the  remainder  of  this  report  concern  steel  shrink- 
fit  assemblies  (E  =  30  x  10^  psi,  U  *  .3)  with  an  assumed  friction  coefficient 
of  0.15. 

Primitive  Model  of  3-D  Shaft  and  Sleeve  Shrinkfit  with  Cyclic  Bending 
Case  1.  The  response  of  the  coarsely  discretized  3-D  shrinkfit  model  to  cyclic 
bending  is  briefly  discussed.  Although  the  predictions  for  a  more  finely  discre¬ 
tized  model  are  more  important,  these  results  demonstrate  that  the  choice  of  fric¬ 
tional  stiffness  is  important  in  cyclic  friction  analysis  by  the  ABAQUS  program. 

The  ABAQUS  shrinkfit  and  bending  idealizations  for  frictional  stiffnesses 
»  2  x  10^  and  6  x  10^  were  subjected  to  several  fully  reversed  bending  load 
cycles  using  the  restart  features  of  ABAQUS.  The  bending  load  amplitude  was  1000 
in-lb.  A  steady-state  condition  was  reached  for  both  stiffnesses  after  only  a 
few  cycles.  The  predicted  sequence  of  slip  events,  normal  stress  vs.  friction 
stress,  friction  stress  vs.  relative  tangential  displacement,  and  relative  tangen¬ 
tial  displacement  vs.  bending  load  are  shown,  at  the  most  crucial  node  pairs, 
for  «  2  x  10^  in  Figures  50-53.  Corresponding  predictions  for  the  »  6  x  10^ 
appear  in  Figures  54-57. 

A  fully  linear  limit  cycle  is  approached  for  the  lower  frictional  stiffness. 

All  slip  eventually  ceases  and  the  assembly  shakes  down  to  a  fully  adhered  state 
within  just  two  cycles.  The  contact  surface  "hardens"  very  quickly  with  this 
particular  combination  of  variables.  The  behavior  of  the  higher  stiffness  system 
is  fundamentally  different  in  that  a  steady-state  slipping  cycle  is  reached  after 
an  initial  transient.  The  contact  surface  does  not  harden  sufficiently  to  prevent 
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a  steady-state  hysteresis  loop  in  which  the  work  done  on  the  system  by  the  bending 
is  dissipated  partially  by  repeated  frictional  slipping.  The  behavior  for  the 
higher  stiffness  is  much  more  interesting  and  probably  more  representative  of  fre- 
ting  conditions.  The  higher  behavior  is  more  fully  investigated  in  a  shrinkfit 
model  with  improved  discretization. 

Irp roved  Model  of  3-D  Shaft  and  Sleeve  Shrinkfit 
with  Cyclic  Bending 

Case  2.  ABAQUS  analysis  of  the  more  finely  discretized  3-D  shrinkfit  model  in 
Figure  43  showed  that  the  major  differences  in  steady-state  slip  predictions 
observed  in  the  coarsely  discretized  model  (fully  adhered  vs.  mixed  adherence/slip 
limit  cycles)  are  due  more  to  finite  element  discretization  than  to  the  frictional 
stiffness.  This  point  is  clarified  in  Figure  58.  This  figure  shows  that  the  fric¬ 
tional  stiffnesses  producing  steady-state,  fully  adhered  and  partially  adhered 
conditions  in  the  coarse  model  (2  x  10^  and  6  x  10^ ,  respectively)  do  not  predict 
fundamentally  different  slip  behaviors  when  the  finite  element  mesh  is  made  finer 
in  the  radial  and  axial  directions.  Clearly  then,  when  the  frictional  stiffness  is 
applied  to  analyze  these  complex  frictional  systems,  simultaneous  mesh  convergence 
and  parametric  frictional  stiffness  studies  are  needed. 

The  steady-state  mixed  adherence/slip  conditions  depend  heavily  on  the  bending 
moment,  as  shown  in  Figure  59.  In  all  cases  frictional  slip  was  principally  con¬ 
fined  to  the  three  percent  of  contact  zone  length  closest  to  the  sleeve  edge; 
however,  the  results  for  1000  in-lb.,  666-2/3  in-lb.  and  333-1/3  in-lb.  maximum 
bending  moment  magnitudes  indicate  that  frictional  slip  becomes  more  severe  and 
widespread  as  bending  moment  increases,  a  result  qualitatively  in  agreement  with 
the  DTNSRDC  tests.  A  rough  quantitative  comparison  of  theory  and  experiment  in 
Table  7  shows  that  the  combination  of  discretization  (36  3-D  elements,  15  contact 
node  pairs),  frictional  stiffness  (K^  -  6  x  10^),  and  frictional  coefficient 
(  p  ■  0.15)  used  in  the  analyses  may  have  produced  pessimistic  slip  predictions. 
This  illustrates  the  not  unsurmountable  difficulty  in  predicting  complex  fric¬ 
tional  behavior  in  mechanical  systems  by  a  method  that  requires  knowledge  of  a 
contact  surface  stiffness  measure  that  only  grossly  represents  many  influential 
tribological  variables. 


Figures  60-62  show  the  relation  of  frictional  shear  stress  to  normal  contact 
stress  at  the  two  most  highly  loaded  node  pairs  for  the  three  bending  load  levels. 
The  stress  behavior  at  these  two  points  approaches  steady-state  loops  that  are 
identical  to  each  other  but  180  deg  out  of  phase.  The  loops  trace  progressively 
smaller  excursions  in  the  normal  force  direction  as  bending  load  decreases.  In  all 
cases,  the  steady  state  is  quickly  reached  after  a  brief  transient  period.  The 
decrease  in  cycle  duration  spent  under  slip  conditions  with  decrease  in  bending 
load  is  clearly  evident. 

The  predicted  relative  tangential  contact  surface  displacements  vs.  frictional 
shear  stress  are  shown  for  1000  in-lb,  666  2/3  in-lb,  and  333  1/3  in-lb  moments  In 
Figures  63-65.  The  displacement  vs.  stress  behavior  at  the  two  most  highly  loaded 
node  pairs  approaches  hysteresis  loops  that  are  identical  in  shape  but  180  deg  out 
of  phase.  The  steady-state  loops  are  quickly  approached  after  brief  transients 
and  enclose  a  larger  "area"  as  bending  moment  increases.  Not  unexpectedly,  the 
frictional  energy  dissipation  represented  by  this  hysteretic  behavior  becomes  much 
more  pronounced  as  bending  load  increases.  The  energy  dissipated  at  one  node 
pair  per  cycle  can  be  determined  easily.  The  shear  stress  axis  can  be  converted  to 
a  force  axis  by  a  simple  multiplicative  factor,  and  the  area  enclosed  by  a  single 
steady-state  hysteresis  loop  can  then  be  calculated.  The  total  energy  dissipated 
on  the  contact  surface  per  cycle  is  found  by  adding  the  steady-state  hysteresis 
loop  areas  for  all  slipping  node  pairs. 

A  side  issue  not  explored  in  this  analysis  is  the  effect  of  shrinkfitting  stage 
analysis  method  on  subsequent  response  to  mechanical  loading.  All  analyses  have 
treated  the  shrinkfitting  as  a  purely  mechanical  loading,  while  in  reality  the  pro¬ 
cess  is  both  mechanical  and  thermal.  (The  sleeve  is  heated  to  high  temperature  and 
is  slipped  onto  the  shaft  while  thermally  expanded.  The  sleeve  shrinks  tightly  onto 
the  shaft  as  it  cools).  This  thermal  process  can  be  simulated  with  ABAQUS  and  will 
certainly  produce  initial  shrinkfitting  and  transient  cycling  response  that  differs 
from  that  presented  here.  It  is  believed,  however,  that  the  thermomechanical  option 
of  shrinkfit  modelling  and  the  purely  mechanical  approach  used  here  will  produce 
similar  steady-state  frictional  behavior. 
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APPLICATION  OF  LAGRANGE  MULTIPLIER  METHOD  TO  CONTACT/GAP  PROBLEMS 

The  mixed  contact/gap  analysis  capability  of  ABAQUS  has  not  been  exercised  in 
this  work.  It  is  currently  being  applied,  in  a  separate  project,  to  analysis  of  the 
interaction  of  controllable  pitch  propeller  blade  components.  The  ABAQUS  results 
obtained  to  date  seem  reasonable  and  are  computed  quite  straightforwardly  as  long 
as  contact  surface  linearization  constraints  are  used.  The  results  of  these  efforts 
will  be  published  in  a  separate  DTNSRDC  technical  memorandum. 

POSSIBLE  IMPROVEMENTS 

IMPROVED  CONTACT  ELEMENTS 

The  aforementioned  experiences  in  contact  analysis  revealed  a  need  for  contact 
restraints  that  are  mathematically  compatible  with  higher-order  finite  elements. 
Although  this  difficulty  was  circumvented  here  by  contact  surface  displacement 
linearizations,  such  a  mathematical  artifice  detracts  from  solution  accuracy.  This 
problem  can  be  solved  once  and  for  all  by  "interface  elements"  based  on  assumed 
contact  pressure  and  contact  surface  displacement  interpolations  that  are  mathemati¬ 
cally  compatible  with  2-D  and  3-D  higher-order  solids. 

A  family  of  such  elements  has  been  derived  for  use  in  conjunction  with  first- 
and  second-order  solids-^*0  These  interface  elements  are  based  on  the  concept  that 
the  contact  surface  is  a  2-D,  geometrically  continuous  "sheet"  of  sorts  over  which 
pressures  and  displacements  also  vary  continuously.  The  conventional  isoparametric 
Interpolation  method  (reference  39,  Chapter  8)  Is  used  in  conjunction  with  trapezoi¬ 
dal  Integration  over  the  surface  in  the  first-order  case  and  Simpson  Rule  integration 
in  the  second-order  case.  In  the  conventional  direct  stiffness  equation  assembly 
approach,  the  interpolated  interface  elements  possess  elemental  contact  force  vec¬ 
tors  that  combine  to  form  the  global  contact  force  vector.  In  contrast,  pointwlse 
contact  node  pair  constraints  possess  individual  contact  force  components.  Unlike 
the  Independent  forces  of  the  pointwlse  contact  element  approach,  the  elemental 
contact  force  vectors  are  inherently  compatible  with  both  first-  and  second-order 
solids.  Interface  element  "strains"  are  defined  as  the  relative  displacements 
between  contact  node  pairs.  These  displacements  are  monitored  for  contact  and  gap¬ 
ping.  The  active  node  pairs  in  contact  are  imposed  via  the  Lagrange  multiplier 
technique  while  inactive  (gapping)  node  pairs  are  determined  by  iteration. 


In  summary,  the  definition  of  the  global  contact  force  vector  is  the  one  major 
difference  between  this  new  Interface  approach  and  the  pointwise  contact  element 
method  used  herein.  Additional  details  on  the  interpolated  interface  theory  and 
results  produced  in  some  simple  test  cases  with  the  ABAQUS  code  are  given  in  refer¬ 
ence  38.  The  interpolated  interface  method  performs  quite  well  in  simple  2-D  and 
3-D  Hertzian  contact  test  problems. 

IMPROVED  FRICTION  MODELLING 

These  experiences  in  friction  analysis  have  shown  that  the  frictional  stiff¬ 
ness  method  is  computationally  useful  and  is  capable  of  producing  reasonable 
results.  It  is  not  always  easy  to  choose  a  proper  frictional  stiffness.  The 
widely  varying  results  produced  by  different  finite  element  discretizations  with 
the  same  frictional  stiffness  suggest  that  frictional  stiffness  is  a  problem- 
dependent  parameter  as  much  as  a  material  property  parameter.  In  addition,  the 
relation  of  this  stiffness  to  tribological  conditions  is  vague  and  has  physical 
meaning  only  in  the  sense  of  the  resistance  of  surface  asperities  in  shear.  An 
effort  has  been  made  to  identify  a  promising  new  approach  for  metallic  friction 
based  on  experimentally  definable  variables  that  better  reflect  the  resistance  of 
contact  surfaces  to  deformation.  One  such  method  is  briefly  described  in  the 
following  paragraphs. 

This  new  approach  to  contact  is  called  "Critical  State  Theory"  (CST)  which 
refers  to  a  mathematical  model  for  soil  mechanics  from  which  it  borrows  some  con¬ 
cepts.  CST  abandons  the  mathematically  convenient  Lagrange  multiplier  artifice  for 
normal  contact  in  favor  of  a  nonlinear  normal  pressure-normal  strain  relation. 

This  relation  produces  experimentally  observed  "hardening"  and  "softening"  behavior 
through  a  pressure-dependent  nonlinear  stiffness.  The  gradual  approach  of  contact 
surface  asperities  during  loading  and  the  gradual  separation  of  asperities  in  un¬ 
loading  is  well  represented  by  the  mathematical  "hardening”  and  "softening"  behavior. 
The  stiffness  in  shear  is  represented  by  a  conventional  shear  modulus  which  utilizes 
the  nonlinear  normal  pressure-normal  strain  relationship  instead  of  the  usual  elas¬ 
tic  modulus.  The  shear  behavior  is  thus  coupled  to  the  normal  behavior  and  reflects 
a  hardening  shear  stiffness  with  increasing  normal  pressure.  This  relation  repre¬ 
sents,  in  a  gross  way,  the  increasing  asperity  resistance  to  shear  loading  that 
occurs  with  Increasing  clamping  pressure.  Since  they  are  valid  in  the  sub-slip 
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range,  where  all  strains  are  recoverable,  the  normal  and  shear  stiffnesses  just 
discussed  are  nonlinear  but  elastic.  Both  of  these  stiffness  formulations  are  thus 
meaningful  in  a  tribological  sense  and  are  explicitly  defined  in  terms  of  experimen¬ 
tally  obtainable  variables. 

The  CST  approach  for  frictional  slip  strongly  resembles  classical  strain¬ 
hardening  flow  theories  of  plasticity.  Some  features  of  plasticity  theories  are 
assumed  in  the  definition  of  slip  displacements;  namely,  strain  rate  decomposition, 
associated  flow,  normality  rule,  and  the  yield  surface  concept.  The  hardening  rule 
is  a  form  of  combined  isotropic  and  kinematic  hardening  in  which  the  size  and  orien¬ 
tation  of  the  "yield"  (in  this  case,  slip)  surface  depends  on  the  current  relative 
tangential  strain.  The  slip  surface  represents  a  critical  intensity  of  contact 
surface  shear  and  normal  forces  at  which  s'lip  initiates  and  continues.  The  extent 
to  which  the  classical  Coulomb  friction  limit  (analagous  to  perfectly  plastic  flow) 
is  reached  depends  on  the  hardening  parameters.  The  "plastic  strain"  components, 
calculated  by  conventional  normality  and  associated  flow  assumptions,  represent 
nonrecoverable  contact  surface  displacements  that  occur  under  slipping  conditions. 
These  "plastic  strains"  are,  in  a  sense,  tribologically  real  because  slip  occurs 
when  contact  asperities  are  grossly  (plastically)  deformed  and  forced  over  one 
another  by  shear  loads.  The  geometric  scale  of  this  predicted  frictional  plasticity 
is  correctly  confined  to  the  contact  surface.  The  CST  frictional  slip  formulation 
is  also  capable  of  producing  the  energy  dissipation  by  hysteresis  observed  in 
cyclically  loaded  frictional  systems. 

QO 

In  summary,  the  CST  concept  defines  contact  surface  behavior  in  meaningful 
tribological  terms  rather  than  relying  on  mathematical  conveniences.  This  theory 
is  one  step  toward  a  rational  description  of  interface  behavior  in  "constitutive 
equation”  terms  and  shows  promise  as  a  basis  for  further  generalization  toward 
certain  kinds  of  lubricated  contacts. 

SUMMARY  AND  CONCLUSIONS 

This  work  has  verified,  extended,  and  improved  finite  element  methods  for  solu¬ 
tion  of  several  classes  of  dry  contact  problems  in  solid  mechanics.  The  class  of 
problems  treated  can  be  categorized  as  follows: 

(1)  Dry  contact,  no  lubrication 

(2)  Static  contact,  no  dynamic  impact  effects 
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(3)  Monotonically  and  quasi-static  cyclically  loaded  contacts 

(4)  Elastic  material  behavior 

(5)  Small  displacements  and  strains 

(6)  Frictionless  contacts  and  contacts  with  Coulomb  friction  under 
both  adherence  and  mixed  slip  adherence  conditions 

(7)  Full  contact  (Mixed  contact  gapping  has  been  treated  in  a  separate 

effort) 

The  contact  analysis  capabilities  of  the  MARC  and  ABAQUS  programs  have  been 
verified  by  comparison  of  elasticity  solutions  to  computed  responses  of  contacts 
with  and  without  friction.  The  Lagrange  multiplier  method  of  MARC  and  ABAQUS  han¬ 
dles  the  normal  contact  problem  quite  satisfactorily.  The  stiffness  method  is  best 
suited  to  the  friction  problem.  Because  the  ABAQUS  program  follows  these  two 
approaches,  it  is  recommended  as  an  analytical  tool  for  unlubricated  contacts. 

The  existing  ABAQUS  capability  has  been  extended  to  modelling  with  higher-order 
elements  and  has  been  applied  to  problems  that  must  be  posed  in  three  dimensions. 

The  current  method  of  modelling  with  higher-order  elements  is  somewhat  deficient 
because  of  limitations  imposed  by  necessary  compatibility  restraints.  The  new  con¬ 
tact  modelling  concept  of  reference  41  could  not  be  evaluated  In  this  work,  but 
preliminary  tests  have  shown  that  It  will  effectively  eliminate  contact  constraint/ 
higher-order  solid  compatibility  difficulties. 

ABAQUS  has  been  applied  to  a  3-D  contact  problem  which  represents  most  features 
of  a  fretting  corrosion  and  fretting  fatigue  test  apparatus.  The  results  obtained 
for  this  study  demonstrate  that  the  complex  contact  and  friction  interaction  occurring 
under  fretting  conditions  can  be  modelled  mathematically.  The  contact  surface  energy 
dissipation  and  the  portion  of  dissipated  energy  leading  to  surface  damage  can  be 
calculated  by  finite  element  methods. 

Despite  Its  limitations,  this  work  has  verified  advanced  methods  for  analyzing 
unlubricated  mechanical  contacts  and  has  extended  the  analysis  capability  to  a  level 
that  is  useful  in  complex  practical  problems. 

TOPICS  FOR  FUTURE  WORK 

This  experience  has  demonstrated  that  there  are  many  "side  issues"  peculiar 
to  certain  dry  contact  problems  which  can  only  be  resolved  after  prolonged  effort. 
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These  side  issues  were  not  thoroughly  treated  in  this  work,  but  they  can  be  con¬ 
sidered  as  open  topics  worthy  of  further  research. 

(1)  If  singular  points  exist  in  the  stress  field  due  to  sharp  corners  or 
cracks,  computed  stresses  will  not  be  realistic  in  the  vicinity  of  the  corner  or 
crack  unless  the  singularity  is  explicitly  modelled.  Such  stress  concentrations 
may  also  cause  material  nonlinearities  in  the  form  of  plasticity.  All  of  this  will 
have  some  effect  on  contact  surface  behavior. 

(2)  The  effect  of  element  size  and  disposition  on  the  convergence  of  displace¬ 
ment  and  stress  predictions  toward  some  "actual”  result  should  always  be  examined 

in  finite  element  treatments  of  nonlinear  problems.  This  issue  could  not  be  fully 
treated  in  the  many  test  and  evaluation  cases  of  this  study.  Although  20-node 
solids  are  monotonically  convergent  In  themselves,  convergence  of  an  entire  solid 
element-contact  element  system  in  the  presence  of  singularities  (if  any),  contact 
surface  constraints  (if  any),  and  plasticity  (if  any)  is  not  guaranteed  and  should 
be  checked. 

(3)  Matters  are  not  simple  if  stresses  very  near  the  contact  surface  are  of 
interest.  Extrapolation  of  integration  point  stresses  proved  useful  in  this  study, 
but  there  Is  no  single  foolproof  way  to  interpret  the  discontinuous  stress  fields 
that  naturally  occur  in  displacement-based  finite  elements.  Alternate  means  of  cal¬ 
culating  surface  stress  should  be  evaluated.  One  useful  alternative  may  be  hybrid 
finite  elements,  which  are  mathematically  constructed  to  predict  continuous  stress 
fields.  Hybrid  elements  also  converge  faster  than  displacement-based  elements  in 
zones  adjacent  to  a  singularity Some  very  recent  contact  analysis  work  in  this 
direction  has  been  accomplished.^ 

A  more  fundamental  problem  involves  unlubricated  interface  modelling.  The  con¬ 
tact  and  friction  behavior  of  interfaces  has  been  simulated  by  necessarily  simple 
mathematical  conveniences  (Lagrange  multipliers  and/or  constraint  equations),  but 
such  models  are  somewhat  artificial.  In  many  cases,  such  approximate  models  can 
capture  the  essentials  of  contact  interface  behavior,  but  this  holds  true  only  in 
cases  where  friction  stresses  are  limited  by  the  Coulomb  theory.  In  particular, 
although  the  stiffness  approach  to  friction  is  computationally  successful,  a  fric¬ 
tional  stiffness  characteristic  that  in  some  way  reflects  surface  roughness  is 
required.  The  computed  results  will  depend  on  the  chosen  stiffness,  and  there  seems 
to  be  no  clear  way  of  correlating  this  number  to  actual  surface  characteristics. 
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There  is  a  need  for  experimentally  definable  contact  surface  friction  models,  ex¬ 
pressed  in  terms  of  finite  element  constitutive  equations,  that  more  directly 
account  for  surface  roughness.  Some  possible  first  steps  in  this  direction  have 
been  discussed. 

The  finite  element  method  is  a  very  powerful  analytical  tool  and  offers  promise 
for  treatment  of  even  more  complex  interface  problems  in  machines  and  machine  com¬ 
ponents.  An  example  of  some  advanced  dry  contact  work  is  a  finite  element  Simula- 

ft 

tion  of  transient  thermoelastic  contact  with  wear  effects.  Lubricated  contacts 
are  another  important  class  of  problem;  much  work  has  been  done  in  developing 
finite  element  methods  to  determine  pressure  distributions  and  load  capacities  for 
’complex  bearings  operating  in  fully  lubricated  conditions For  highly  loaded 
bearings  or  certain  bearing  pad  materials,  the  bearing  pad  deformations  are  as 
important  an  effect  as  the  lubricant  behavior.  Finite  element  methods  have  been 
used  to  calculate  pressures,  film  thickness  distributions,  and  maximum  load  capaci¬ 
ties  of  bearings  under  elastohydrodynamic  conditions. ^  ^1  Finally,  finite  ele¬ 
ments  have  been  applied  to  lubricated  contacts  in  which  thermal  effects  are 
72 

important.  It  appears  that  well-defined  numerical  methods  exist  for  analysis 
of  geometrically  complex,  fully  lubricated  contacts  for  a  variety  of  conditions. 
Although  lubricated  contacts  are  very  important  in  naval  machinery  applications, 
such  problems  require  analytical  approaches  which  in  no  way  resemble  those  con¬ 
sidered  herein  for  dry  contact.  Finite  element  treatment  of  such  problems  would 
require  completely  new  efforts. 

Finally,  the  problem  of  partially  lubricated  contacts  (boundary  lubrication) 
is  much  more  complex  than  dry  or  fully  lubricated  cases.  None  of  the  dry  contact 
methods  considered  here  or  the  lubricated  contact  approaches  listed  above  are  appro¬ 
priate.  Application  of  finite  elements  to  boundary  lubrication  first  requires  a 
statement  of  the  problem  in  fundamental  mathematical  terms.  Such  mathematical  state¬ 
ments  must  reflect  observed  relationships  among  important  variables,  which  can  only 
be  defined  through  experimental  effort.  Even  so,  finite  elements  may  not  be  the 
most  effective  technique  for  analyzing  boundary  lubrication  problems. 
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Contact  Forces  Resulting  from  Uniform  Pressure  on  One  Face  of  Rectangular  3-D  Solids 


Figure  5  -  Elastic  Shaft/Rigid  Sleeve  Contact 


Figure  6  -  Finite  Element  Representation  of  Elastic  Shaft/Rigid  Sleeve  Contact 
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Figure  11  -  Elastic  Shaft/Rigid  Sleeve,  Axial  Shaft  Deflections  at  Inner  and  Outer  Radii 


Figure  13  -  Elastic  Shaft/Rigid  Sleeve,  Hoop  Stresses  Near  Contact 
Surface  and  in  Shaft  Interior 


Figure  13a  -  Hoop  Stress  Near  Contact  Surface 


Figure  13  (Continued) 
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Figure  13c  -  Stress  at  7/10  Outer  Radius 
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Figure  14  -  Elastic  Shaft/Rigid  Sleeve,  Axial  Stress  Near  Shaft  Surface 


510 f-  508  p*i 


Figure  16  -  Elastic  Shaft/Rigid  Sleeve,  Contact  Element  Pressures 
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Figure  17  -  Elastic  Shaf t/Elastic  Sleeve  Contact 


Figure  18  -  Finite  Element  Representation  of  Elastic  Shaft/Elastic  Sleeve  Contact 
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Figure  20  -  Typical  Stress  Prediction  of  Fully  Integrated  Higher  Order  Solids  Near  a 
Steep  Stress  Gradient 
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Elastic  Shaft/Elastic  Sleeve,  Axial  Displacements  at  Contact  Surface 


CONTACT  ZONE 
HALF  LENGTH 


Figure  25  -  Elastic  Shaft/Elastic  Sleeve,  Radial  Stresses  Extrapolated  to  Contact  Surface 
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Elastic  Shaft/Elastic  Sleeve,  Axial  Stress  Extrapolated  to  Contact  Surface 


ire  28  -  Elastic  Shaft /Elastic  Sleeve,  Contact  Element  Pressure  Compared  to  Extrapolated 
Radial  Stress 
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Figure  30  -  Elastic  Shaft/Elastic  Sleeve,  Comparison  of  Radial  Deflection  Predictions  for 
Axisymmetrlc  and  3-D  Models 
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ELEMENT  1  ELEMENT  2  CORNER  NODE  CONTACT  PRESSURE 


Figure  33  -  Contact  Pressures  Calculated  from  Contact  Forces 


Figure  34  -  Elastic  Shaft/Elastic  Sleeve.  Comparison  of  Contact  Element  Pressure  for  Axisymmetric 
and  3-D  Models 
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NOTE:  ALL  ELEMENTS  IN  FULL  ADHERENCE 
IX  *  0.15  ASSUMED 


Figure  39  -  Coarse  Finite  Element  Mesh  for  3-D  Elastic  Shaft/Elastic  Sleeve 
Shrinkfit  with  Bending  and  Friction 
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FRETTING  TEST  SETUP 


SHAFT 

COLLAR  / 


BENDING  MOMENT 


EQUIVALENT  STRUCTURAL  MODEL 

Figure  40  -  Experimental  Fretting  Corrosion  Test  Machine 
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SYM  SYM  SYM  SYM 


1  RADIAL  SHRINKFIT 

INTERFERENCE  -  0.0005  in.-t> 

NO  BENDING 
NO  SLIP  ANYWHERE 


SHRINKFIT  PLUS 
M  *  333  1/3  m. -lb 

THIS  NODE  PAIR  SLIPPED 
IN  INTERVAL 

222  M  333  in. -lb 


SHRINKFIT  PLUS 
M  -  666  2/3  in.-lb 


THIS  NODE  PAIR 
SLIPPED  IN  INTERVAL 
333  M  •  444  in.-lb 


SHRINKFIT  PLUS 
M  =  1000  in.-lb 


THIS  NODE  PAIR  SLIPPED 
IN  INTERVAL 
889  <  M  <  1000  in.-lb 


Figure  41  -  3-D  Elastic  Shaf t/Elastic  Sleeve  Shrinkfit  with  Bending  and  Friction 
Frictional  Slip  Development  in  Plane  of  Bending  with  Increasing 
Bending  Load 
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Figure  45  -  Coarsely  Discretized  3-D  Shaft/Sleeve  Shrinkfit  with  Bending  and  Friction,  ABAQUS 
Frictional  Slip  Predictions  for  Various  Frictional  Stiffnesses 
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Figure  54  -  Coarse  3-D  Shrinkfit  with  Friction  and  Bending,  Cyclic  Slip  Events 
at  Crucial  Locations  for  ■  6  x  10^ 
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Figure  55a  -  Response  of  Node  Pair  1  (Fig.  50) 
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Figure  55b  -  Response  of  Node  Pair  2  (Fig.  50) 

Figure  55  -  Coarse  3-D  Shrinkfit  with  Friction  and  Bending,  Cyclic  Shear  Versus 
Normal  Stress  at  Crucial  Locations  for  ■  6  x  10^ 
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Figure  56  -  Coarse  3-D  Shrinkfit  with  Friction  and  Bending  Cyclic  Frictional  Shear 
Versus  Relative  Tangential  Displacements  at  Crucial  Locations  for 
K  «  6  x  105 
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1.3  1.45 

Figure  59a  -  Figure  59b  -  Figure  59c  - 

Max  Moment  Max  Moment  Max  Moment 
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Figure  59  -  Fine  3-D  Shrinkflt  with  Friction  and  Bending, Comparison  of  Slip 
Predictions  for  Three  Levels  of  Maximum  Bending  Moment  for 
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Figure  60b  -  Response  of  Node  Pair  2  (Fig.  58) 


Figure  60  -  Fine  3-D  Shrinkfit  with  Friction  and  Bending,  Cyclic  Shear  Versus 
Normal  Stress  at  Crucial  Locations,  *  6  x  10-*,  Maximum  Bending 


Load  ■  1000  Inch-Pounds 
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Figure  63  -  Fine  3-D  Shrinkfit  with  Friction  and  Bending,  Cyclic  Frictional  Shear 
versus  Relative  Tangential  Displacements  at  Crucial  Locations,  “ 

6  x  105,  Maximum  Bending  Load  *  1000  Inch-Pounds 
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Figure  64a  -  Response  of  Node  Pair  1  (Fig.  58) 
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Figure  64b  -  Response  of  Node  Pair  2  (Fig.  58) 

Figure  64  -  Fine  3-D  Shrinkfit  with  Friction  and  Bending,  Cyclic  Frictional  Shear 
versus  Relative  Tangential  Displacements  at  Crucial  Locations,  Kf  » 

6  x  10^,  Maximum  Bending  Load  *  666  2/3  Inch-Pounds 


Figure  65a  -  Response  of  Node  Pair  1  (Fig.  58) 


Figure  65b  -  Response  of  Node  Pair  2  (Fig.  58) 


Figure  65  -  Fine  3-D  Shrinkfit  with  Friction  and  Bending,  Cyclic  Frictional  Shear 
Versus  Relative  Tangential  Displacements  at  Crucial  Locations,  Kf  - 
6  x  105,  Maximum  Bending  Load  -  333  1/3  Inch-Pounds 
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TABLE  1  -  RESULTS  OF  MARC  CALCULATIONS;  ATTEMPTS  TO  SIMULATE 
CONWAY'S  RIGID  PUNCH/ELASTIC  SLAB  FRICTIONAL  CONTACT  PROBLEM 


Punch  Half  Width:  Slab  Depth  ■  1:4 

Approximate  % 

of  Adherence 

Validity 

Friction 

(from  Ref.  49,  see 

MARC  Result 

of  MARC 

Coefficient 

Figure  38) 

Solution 

0.60 

99 

Full  Adherence  -  Converged 

Valid 

0.50 

97 

Full  Adherence  -  Slow  Convergence 

NV* 

0.40 

90 

Full  Slip  -  Slow  Convergence 

NV 

0.30 

78 

Full  Slip  -  No  Convergence 

NV 

0.20 

55 

Partial  Slip  -  No  Convergence 

NV 

0.10 

14 

Full  Adherence  -  Slow  Convergence 

NV 

Punch  Half  Width:  Slab  Depth  =  1:2 

0.60 

99 

Full  Adherence  -  Converged 

Valid 

0.55 

98 

Full  Adherence  -  Converged 

Valid 

0.50 

97 

Full  Adherence  -  Slow  Convergence 

NV 

0.45 

94 

Full  Adherence  -  Slow  Convergence 

NV 

0.40 

90 

Full  Slip  -  Slow  Convergence 

NV 

0.35 

85 

Full  Slip  -  Slow  Convergence 

NV 

0.30 

78 

Full  Slip  -  No  Convergence 

NV 

0.25 

68 

Full  Slip  -  No  Convergence 

NV 

0.20 

55 

Partial  Slip  -  No  Convergence 

NV 

0.15 

33 

Partial  Slip  -  No  Convergence 

NV 

0.10 

14 

Full  Adherence  -  Slow  Convergence 

NV 

0.05 

4 

Full  Adherence  -  Slow  Convergence 

NV 

*NV  =  Not  Valid 
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TABLE  2  -  RESULTS  OF  MARC  CALCULATIONS;  CONTACT  SURFACE  BEHAVIOR  OF  3-D 
ELASTIC  SHAFT/ELASTIC  SLEEVE  SHRINKFIT  WITH  BENDING  AND  FRICTION 


Ratio  Fip/Fjj  (Tangential  Contact  Force/Normal  Contact  Force)  Brackets 

Indicate  Slipping  Condition 


Node  Pair 
Number* 


Shrinkfit  and 
Bending 

Shrinkfit  and 
333  1/3  in-lb 

Shrinkfit  and 
666  2/3  in-lb 

Shrinkfit  and 
1000  in-lb 

0 

0 

0 

0 

0.0023 

0.0059 

0.0879 

0.1275 

0.0184 

0.1188 

[0.1491] 

[0.1501] 

0 

0 

0 

0 

0.0017 

0.0403 

0.0455 

0.0523 

0.0174 

0.0549 

0.0972 

0.1433 

0 

0 

0 

0 

0.0023 

0.0392 

0.1192 

[0.1517] 

0.0184 

[0.1645] 

[0.1480] 

[0.1500] 

%  Shrinkfit  Restraint  (Percentage  of  Prescribed  Shrinkfit 

Restrained  by  Shaft) 


86 

85 

85 

84 

85 

89 

94 

98 

87 

99 

112 

125 

86 

86 

86 

86 

85 

85 

85 

85 

87 

87 

87 

87 

86 

65 

65 

64 

85 

60 

55 

49 

87 

54 

40 

25 

♦Location  of  Node  Pairs. 


PI 

m 


1 


TABLE  3  -  ABAQUS  ANALYSES  OF  RIGID  PUNCH/ELASTIC  SLAB  CONTACT 
WITH  FRICTION,  n  -  0.35 


Ratio  F'p/F^j  (Tangential  Contact  Force/Normal  Contact 
Brackets  Indicate  Slipping  Condition 

Force) 

Distance  from 
Center  Plane 

Frictional  Stiffness  (K^) 

of  Symmetry 
(in. ) 

3xl06 

lxlO7 

l.lxlO7 

1 . 36xl07 

1 . 37xl07 

5xl07 

0.0 

0 

0 

0 

0 

0.250 

0.017 

0.019 

0.019 

0.019 

0.375 

0.040 

0.044 

0.044 

0.045 

0.500 

0.625 

0.054 

0.070 

0.063 

0.087 

0.064 

0.088 

0.065 

0.092 

Noncon- 

vergent 

Noncon- 

vergent 

0.750 

0.114 

0.142 

0.144 

0.150 

Solution 

Solution 

0.8125 

0.179 

0.226 

0.231 

0.242 

0.875 

0.197 

0.276 

0.285 

0.311 

0.9375 

0.212 

0.341 

[0.350] 

[0.350] 

1.0 

0.196 

[0.350] 

[0.350] 

[0.350] 

no 
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TABLE  4  -  ABAQUS  ANALYSES  OF  RIGID  PUNCH/ELASTIC  SLAB 
CONTACT  WITH  FRICTION,  |J  =  0.80 


Ratio  F^/F^  (Tangential  Contact  Force/Normal  Contact 
Force)  Brackets  Indicate  Slipping  Condition 


Distance  From 
Centerplane  of 
Symmetry  (in.) 

Frictional  Stiffness  (Kf) 

l.lxlO7 

1. 14xl07 

1.15xl07 

0.0 

0 

0 

0.250 

0.018 

0.018 

0.375 

0.041 

0.041 

0.500 

0.060 

0.059 

Noncon- 

Noncon- 

0.625 

0.083 

0.083 

vergent 

vergent 

0.750 

0.135 

0.134 

Solution 

Solution 

0.8125 

0.214 

0.214 

0.875 

0.262 

0.263 

0.9375 

0.319 

0.320 

1.0 

0.424 

0.431 

TABLE  5  -  ABAQUS  ANALYSES  OF  RIGID  PUNCH/ELASTIC  SLAB 
CONTACT  WITH  FRICTION,  FRICTIONAL  STIFFNESS  »  I.lxlO7 


Ratio  F^/Fj,j  (Tangential  Contact  Force/Normal  Contact  Force) 
Brackets  Indicate  Slipping  Condition 


Friction  Coefficient  (^  ) 


Distance  From 
Center  Plane 
of  Symmetry 
(in.) 


0.0 
0.250 
.375 
.500 
0.625 
0.750 
0.8125 
0.875 
0.9375 
1.0 


0.10 

0.15 

0.20 

0.35 

0.50 

0 

0 

0 

0 

0 

0.053 

0.036 

0.029 

0.019 

0.018 

[0.100] 

0.132 

0.067 

0.044 

0.041 

[0.100] 

0.128 

0.099 

0.064 

0.060 

[0.100] 

[0.150] 

0.143 

0.088 

0.083 

[0.100] 

[0.15C] 

[0.200] 

0.144 

0.135 

[0.100] 

[0.150] 

[0.200] 

0.231 

0.214 

[0.100] 

[0.150] 

[0.200] 

0.285 

0.262 

[0.100] 

[0.150] 

[0.200] 

[0.350] 

0.319 

[0.100] 

[0.150] 

[0.200] 

[0.350] 

0.424 

0 

0.018 
.041 
.060 
.083 
.135 
.214 


0.262 

0.319 

0.424 


TABLE  6  -  COMPARISON  OF  ABAQUS  AND  MARC  ANALYSES  OF  COARSELY  DISCRETIZED 
3-D  SHAFT/SLEEVE  SHRINKFIT  WITH  BENDING  AND  FRICTION 


Ratio  Ft/F^  (Tangential  Contact  Force/Normal  Contact  Force) 
Axisymmetric  Shrinkfit  Only 

Node  Pair 

Number* 

MARC 

ABAQUS 

Frictional  Stiffness 

(Kf ) 

2x10* 

mm 

1 

0 

0 

0 

0 

0 

0 

2 

0.0023 

0.0019 

0.0038 

0.0054 

0.0071 

0.0087 

3 

0.0184 

0.0081 

0.0161 

0.0230 

0.0303 

0.0372 

4 

0 

0 

0 

0 

0 

0 

5 

0.0017 

0.0009 

0.0019 

0.0027 

0.0036 

0.0044 

6 

0.0174 

0.0040 

0.0081 

0.0177 

0.0156 

0.0192 

7 

0 

0 

0 

0 

0 

0 

8 

0.0023 

0.0019 

0.0038 

0.0054 

0.0071 

0.0078 

9 

0.0184 

0.0081 

0.0161 

0.0230 

0.0303 

0.0372 

Ratio  Ft/Fn,  Maximum  Bending  Load  1000  in-lb 

Brackets  Indicate  Slipping  Condition 

1 

0 

0 

0 

0 

0 

0 

2 

0.1275 

0.0010 

0.0021 

0.0031 

0.0041 

0.0051 

3 

[0.1500] 

0.0022 

0.0043 

0.0063 

0.0083 

0.0103 

4 

0 

0 

0 

0 

0 

0 

5 

0.0523 

0.0009 

0.0019 

0.0028 

0.0036 

0.0045 

6 

0.1433 

0.0041 

0.0082 

0.0121 

0.0160 

0.0198 

7 

0 

0 

0 

0 

0 

0 

8 

[0.1517] 

0.0026 

0.0052 

0.0076 

0.0100 

0.0122 

9 

[0.1500] 

0.0157 

0.0310 

0.0456 

0.0597 

0.0732 

*Location  of 

Node  Pairs 

as  Shown  in  Table  2. 
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TABLE  7  -  FINE  3-D  SHRINKFIT  WITH  FRICTION  AND  BENDING, 
COMPARISON  OF  EXPERIMENTAL  FRETTING  DAMAGE  EXTENT  WITH 
CALCULATED  SLIP  ZONE  EXTENT 


Maximum 

Bending 

Moment 

Extent  of  Damage  Observed  in 
DTNSRDC  Rotating  Bending 
Experiment 

Extend  of  Significant 
Slip  Calculated  in 
•  ABAQUS  Cyclic  Bend¬ 
ing  Analysis 

(in. /lb) 

Average 
(in. ) 

Maximum 
(in. ) 

(in.) 

333  1/3 

None 

None 

At  Edge  Only 

666  2/3 

0.007 

0.010 

0.05 

1000 

0.009-0.011 

0.018 

0.05 
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APPENDIX 

EDGE  EFFECTS  IN  CONTACT  PROBLEMS 


An  edge  effect,  in  the  form  of  a  wiggly  radial  streas  distribution,  appears 
in  the  extrapolated  MARC  stress  predictions  for  the  elastic  shaft/elastic  sleeve 
shrinkflt  problem.  This  same  phenomenon  was  observed  in  a  corresponding  elastic- 
ity  solution. ^  The  edge  effect  was  found  to  be  more  pronounced  as  the  sleeve 
thickness  become  thinner  in  relation  to  shaft  radius.  The  thin  sleeve  is  modelled 
as  a  thin  shell  in  Hill's  work,  and  the  Influence  of  hoop  restraint  allows  the 
use  of  a  mathematical  analogy  to  a  beam  on  an  elastic  foundation  (reference  73, 
pp.  30-33).  This  analogy  gives  rise  to  the  edge  effect. 

This  edge  effect  can  be  explained  even  for  cases  where  the  sleeve  is  not  thin. 

If  the  shaft  and  sleeve  surfaces  are  imagined  as  two  thin  shells,  then  the  radial, 
axial  and  shear  stiffnesses  of  the  components'  interiors  may  be  imagined  as  arrays 
of  attached  discrete  springs-  The  hoop  restraint  of  both  thick  shells  can,  according 
to  the  analogy,  be  represented  by  a  continuous  elastic  foundation.  Both  shaft  and 
sleeve  are  now  beams  on  elastic  foundations,  and  the  total  effects  of  all  spring 
and  foundation  restraint  can  be  schematically  pictured  as  two  arrays  of  springs 
(Figure  70a). 

This  conceptual  assembly  of  beams  on  foundations  is  now  subject  to  the  displac¬ 
ement  condition  that  corresponds  to  the  specified  shrinkflt  interference  (the 
distance  the  sleeve  would  radially  contract  if  unconstrained).  This  results  in  a 
physically  impossible  deformation  as  shown  in  Figure  66a.  Displacement  continuity 
can  be  restored  by  a  concentrated  force  and  moment  at  the  sleeve  edge,  shown  in 
Figure  66b.  The  shear  force  and  bending  moment  patterns  in  beam/foundation  assem¬ 
blies  take  on  the  character  of  damped  waves  when  acted  upon  by  concentrated  forces 
and  moments  (reference  73,  Page  14).  These  disturbances  gradually  die  out  away 
from  the  loaded  point.  This  fact  explains,  in  part,  why  the  stress  distributions 
in  the  3-D  shaft/sleeve  assembly  are  wiggly  near  the  contact  zone  edge.  It  is  a 
natural  consequence  of  the  localized  forces  needed  to  maintain  displacement  con¬ 
tinuity  at  this  edge. 

The  finite  element  procedure  automatically  accounts  for  such  forces  so  that 
both  continuity  and  equilibrium  are  maintained.  The  predicted  displacements  reflect 
some  actual  restrained  shrinkage  compatible  with  both  the  local  and  overall  stiff¬ 
nesses  of  both  shaft  and  sleeve.  Interior  radial  support  of  shaft  and  sleeve 
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augments  the  hoop  restraint  in  the  fictitious  foundation  quite  directly,  but  the 
effects  of  axial  and  shear  restraint  on  the  beam-on-foundation  wiggles  are  difficult 
to  assess.  The  beam-on-foundation  analogy  is  not  perfect,  but  it  at  least  provides 
a  sensible  explanation  for  the  predicted  wiggly  stress  distributions. 
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